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Gq 
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9 

Yt 
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Sf 


St 

Xt 
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1.  INTRODUCTION 


This  report  presents  the  formulation  and  computational  implementation  of  the  LaRC 
ply-based  failure  criteria  for  laminated  composite  materials,  and  the  implementation 
of  a  continuum  damage  model  that  is  able  to  predict  the  onset  and  propagation  of 
ply  failure  mechanisms  as  well  as  the  hnal  failure  load  of  composite  structures. 

The  LaRC  failure  criteria  is  implemented  by  means  of  a  UVARM  ABAQUS®  subrou¬ 
tine,  whereas  the  continuum  damage  model  is  implemented  in  both  a  UMAT  subroutine 
for  implicit  analysis  and  in  a  VUMAT  subroutine  for  explicit  analysis. 

This  report  explains  in  detail  the  dehnition  of  the  required  material  properties  and 
initial  conditions  in  the  ABAQUS®  input  hie.  In  addition,  examples  of  application 
of  the  subroutines  developed  are  presented  at  the  end  of  each  section  of  the  report. 

The  papers  that  were  published  in  the  context  of  the  project  activities  are  pre¬ 
sented  in  Annexes  B  and  C. 
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2.  UVARM  SUBROUTINE 


2.1  Overview 


This  section  presents  a  reference  manual  for  the  use  of  an  ABAQUS®  user  subroutine 
UVARM  [1]  with  the  implementation  of  the  plane-stress  LaRC  failure  criteria  [2]-[5].  The 
basic  equations  of  the  LaRC  failure  criteria  and  the  details  of  its  implementation  are 
presented.  An  example  of  the  use  of  the  subroutine  in  the  simulation  of  an  open-hole 
carbon-epoxy  quasi-isotropic  laminate  loaded  in  tension  is  presented. 


2.2  In-situ  strengths 


The  in-situ  effect,  originally  detected  in  Parvizi’s  tensile  tests  of  cross-ply  glass  hber 
reinforced  plastics  [6] ,  is  characterized  by  higher  transverse  tensile  and  shear  strengths 
of  a  ply  when  it  is  constrained  by  plies  with  different  hber  orientations  in  a  laminate, 
compared  with  the  strength  of  the  same  ply  in  an  unidirectional  laminate.  The  in- 
situ  strength  depends  on  the  number  of  plies  clustered  together,  and  on  the  hber 
orientation  of  the  constraining  plies.  The  model  proposed  for  the  calculation  of  the 
in-situ  strengths  uses  the  simplifying  assumption  that  the  hber  orientation  of  the 
constraining  layers  does  not  ahect  the  value  of  the  in-situ  strength. 

The  in-situ  ehect  is  illustrated  in  Figure  2-1,  where  the  relation  between  the  in- 
situ  transverse  tensile  strength  experimentally  measured  and  the  total  number  of  90° 
plies  clustered  together  (2n)  is  represented. 

It  is  clear  that  accurate  in-situ  strengths  are  necessary  for  any  stress-based  failure 
criterion  for  matrix  cracking  in  constrained  plies.  Therefore,  the  user  should  calculate 
the  in-situ  strengths  that  are  required  for  the  LaRC  failure  criteria  [2]- [5].  The  in-situ 
strengths  are  an  input  for  the  ABAQUS®  UVARM  subroutine. 

The  closed-form  solutions  to  predict  the  in-situ  strengths  previously  proposed  can 
be  used  [2].  The  tensile  transverse  in-situ  strengths  are  given  by  [2]: 


/  8G 

For  a  thin  embedded  ply:  =  -i  /  — : — 

yvrtA^^ 

(2,1) 

/  Gr2+ 

For  a  thin  outer  ply:  Y-p  =  1.79w 

y  T^tA.22 

(2,2) 

For  a  thick  ply:  Yp  =  1.12v^Yp^ 

(2.3) 
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Figure  2-1  In-situ  effect  in  laminated  composites  (after  Dvorak  [7]). 


where  is  the  tensile  transverse  strength  measured  in  an  unidirectional  test  spec¬ 
imen,  t  is  the  ply  thickness,  G2+  is  the  mode  I  fracture  toughness,  and  A22  is  defined 
as: 


A»,  =  2(  A-ta 


I 

E2  El 

The  in-situ  shear  strengths  are  obtained  as  [2]: 


(2.4) 


^L  = 


3/3Gi2 


(2.5) 


where  (3  is  the  shear  response  factor,  and  the  parameter  0  is  defined  according  to  the 
configuration  of  the  ply: 


For  a  thick  ply: 
For  a  thin  ply: 
For  an  outer  ply: 


12  CT' 

Gi2 

48G6 

nt 

24G6 

nt 


uda4 


+ 18/3  (5?“) 


(2.6) 


where  S'l  °  is  the  shear  strength  measured  in  an  unidirectional  test  specimen,  and  Ge 
is  the  mode  II  fracture  toughness. 
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2.3  Failure  criteria 


2.3.1  Transverse  fracture 

Tension 

The  LaRC  criterion  to  predict  failure  under  transverse  tension  ((T22  >  0)  and  in-plane 
shear  is  defined  as: 


/,  \^22,  (  <722 


(1-9) 


a. 


(m) 

22 


Y. 


+  9 


a. 


(m) ' 
22 


Y. 


+ 


(m) 

^12 


2 


1  <  0, 


(Til  <  0,  |(Tii|  <  Xc/2 


(2.7) 


where  g  = 

The  stresses  (jf™^  are  calculated  in  a  frame  aligned  with  the  hber  direction  accord¬ 
ing  to  the  following  expressions: 


(tTi)  9  ,  .9  I  o  I  I  • 

<7ii  =  <7ii  COS  g)  -|-  (T22Sin  99  -|-  2  |(Ti2|  sini/icos^? 
<72™^  =  (Til  sin^  g)  -|-  (T22  cos^  (^  —  2  |(Ti2|  sin  (^  cos  tp 

^(m)  _  gjj^  (.Qg  ^  gjj^  (.Qg 

+  I  <712!  (cos^  (f  —  sin^  if) 
where  the  misalignment  angle  (p  is  dehned  as: 


(2.8) 


<712!  +  (Gi2  —  Xc) 
Gi2 


tan  ^ 


+  <7ll  —  (722 

l-^l-4tu(|;) 

2-^7 


with  CT  =  -I-  7]^. 


(2.9) 


(2.10) 
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Compression 

The  failure  criteria  used  to  predict  fracture  under  transverse  compression  ((T22  <  0) 
and  in-plane  shear  is  dehned  as: 


-  1  <  0,aii  >  -Yc 

(2.11) 

-  1  <  0,aii  <  -Yc 

(2.12) 

The  effective  shear  stresses  in  the  fracture  plane  are  dehned  as: 


=  {\T^\+rf  CFn  cos  0)  (2.13) 

'^e  =  {\r^\+'n^CFnSm.e)  (2.14) 


with  9  =  tan  ^  (^g-JsinL) •  McAuley  operator  dehned  as  {x)  :=  +  |a;|). 

The  components  of  the  stress  tensor  on  the  fracture  plane  are  given  by: 

{O'n  =  <722  COS^  a 

=  — (J22  sin  a  cos  a  (2-15) 

=  <7i2  cos  a 

The  terms  and  are  calculated  from  equations  (2.13)-(2.14)  using  the 
relevant  components  of  the  stress  tensor  established  in  a  frame  representing  the  hbre 
misalignment.  The  fracture  plane  is  dehned  by  the  angle  a.  The  determination  of  a 
is  performed  numerically  maximizing  equations  (2.11)  and  (2.12). 

The  coefficients  of  transverse  and  longitudinal  inhuence,  77^  and  respectively, 
can  be  obtained  as: 


V 


T 


V 


L 


-1 

tan  2q;o 

Sl  cos2q;o 

Yc  cos^  ao 


(2.16) 

(2.17) 


where  ao  is  the  fracture  angle  under  pure  transverse  compression  (ao  ~  53°). 

In  the  absence  of  test  data  the  transverse  shear  strength  can  be  estimated  as: 


St  =  Yc  cosao 


^sin  ao  + 


cosao 

tan2ao 


(2.18) 
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2.3.2  Longitudinal  failure 


Tension 

The  failure  criterion  used  to  predict  fiber  fracture  under  longitudinal  tension  (an  >  0) 
is  defined  as: 

—  -  1  <  0  (2.19) 


Compression 

The  failure  criterion  used  to  predict  hber  fracture  under  longitudinal  compression 
(an  <  0)  and  in-plane  shear  (fiber  kinking)  is  given  as: 


i<^l? 


-  1  <  0,  4?'  <  0 


(1-9) 


a. 


(m) 

22 


Y. 


+  9 


a. 


(m) ' 

22 


Yn 


+ 


a 


(m) ' 
12 


-  1  <  0, 


a. 


(m) 

22 


>0,|an|  >Xc/2 


(2.20) 


Based  on  the  equations  outlined  above,  the  subroutine  calculates  the  following 
parameters  that  can  be  used  for  post-processing: 


Table  2.1  Parameters  calculated  by  the  subroutine. 


UVARM(l) 

Failure  index  for  transverse  tensile  failure 

UVARM(2) 

Failure  index  for  transverse  compressive  failure 

UVARM(3) 

Failure  index  for  longitudinal  tensile  failure 

UVARM(4) 

Failure  index  for  longitudinal  compressive  failure 

2.4  Input  into  Abaqus  standard 


The  user  must  create  a  hie  with  the  name  j  obname .  mt  where  the  material  properties 
are  dehned.  The  hie  must  have  the  same  name  as  the  .  inp  hie  that  dehnes  the  model. 
The  hie  j  obname  .mt  must  be  placed  in  the  same  directory  where  the  ABAQUS®  input 
hie  is  located. 

The  format  of  the  hie  j  obname  .mt  is  the  shown  in  Table  3.1. 

The  symbol  **  means  that  the  corresponding  lines  can  be  used  to  write  comments. 
The  name  of  the  material  (line  4)  must  be  the  same  as  the  one  given  in  the  ABAQUS® 
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Table  2.2  Format  required  for  the  jobname.mt  file. 


Line 

Column 

1 

2 

3 

4 

5 

6 

7 

8 

1 

** 

2 

3 

3 

** 

4 

Material 

5 

** 

6 

El 

L/2 

Ez 

V21 

V32 

7 

** 

8 

Gi2 

^23 

Gai 

Xt 

Xc 

Yt 

Yc 

9 

** 

10 

(To 

/9 

9 

Sl 

jobname.inp  file  and  must  be  written  in  CAPITALS.  Lines  3  to  10  can  be  repeated 
for  the  definition  of  other  materials.  The  following  is  an  example  of  a  file  with  the 
dehnition  of  three  materials. 
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**  LaRC03  failure  criteria:  use  3  for  LaRC03  and  4  for  LaRC04 
3 

**  MAT.  #1:  IM7-8552-thin:  thin  embedded  ply 
IM7-8552-THIN 

**  El  ,  E2  ,  E3  ,  nu21  ,  nu31  ,  nu32 

171420.,  9080.,  9080.,  0.017,  0.017,  0.4 

**  G12  ,  G23  ,  G31  ,  XT  ,  XC  ,  YT  ,  YC  ,  SL 

5290.,  3242.9,  5290.,  2323.5,  1200.1,  160.2,  199.8,  92.3 
**  alphao,  beta,  g  ,  SL_IS 

53.  ,0.  ,  0.5,  130.2 

**  MAT.  #2:  IM7-8552-thin:  thin  outer  ply 
IM7-8552-THIN-0UTER 

**  El  ,  E2  ,  E3  ,  nu21  ,  nu31  ,  nu32 

171420.,  9080.,  9080.,  0.017,  0.017,  0.4 

**  G12  ,  G23  ,  G31  ,  XT  ,  XC  ,  YT  ,  YC  ,  SL 

5290.,  3242.9,  5290.,  2323.5,  1200.1,  101.4,  199.8,  92.3 
**  alphao,  beta,  g  ,  SL_IS 

53.  ,  0.  ,  0.5,  107. 

**  MAT.  #3:  IM7-8552-thin-2t :  embedded  ply  with  t=2*ply  thickness 
IM7-8552-THIN-2T 

**  El  ,  E2  ,  E3  ,  nu21  ,  nu31  ,  nu32 

171420.,  9080.,  9080.,  0.017,  0.017,  0.4 

**  G12  ,  G23  ,  G31  ,  XT  ,  XC  ,  YT  ,  YC  ,  SL 

5290.,  3242.9,  5290.,  2323.5,  1200.1,  113.3,  199.8,  92.3 
**  alphao,  beta,  g  ,  SL_IS 

53.  ,  0.  ,  0.5,  106.9 

The  user  must  define  consistent  material  properties  in  the  jobname.inp  file,  and 
define  four  user  output  variables  following  the  example  shown  below: 

^MATERIAL ,  NAME=IM7-8552-thin-outer 
^ELASTIC,  TYPE=LAMINA 

171420.,  9080.,  0.32,  5290.,  5290.,  3242.9 
*USER  OUTPUT  VARIABLES 
4, 


2.5  Example 

An  ABAQUS®  model  with  an  example  of  the  use  of  the  UVARM  subroutine  in 
the  prediction  of  first  ply  failure  of  a  quasi-isotropic  Hexcel  IM7-8552  [90/0/  ±  45] 3^ 
CFRP  laminate  with  a  central  hole  and  loaded  in  tension  can  be  downloaded  from 
the  following  URL: 

WWW . f e . up . pt/~pcamanho/ oht3_03 . inp 
www.fe.up.pt/~pcamanho/model_oht3 . inp 
WWW . f e . up . pt/~pcamanho/ oht3_03 . mt 
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The  model  uses  ABAQUS®  S4  shell  elements.  The  specimen  is  3mm  thick,  150mm 
long,  36mm  wide  and  has  a  central  hole  with  a  diameter  of  6mm. 

The  material  properties  used  are  shown  in  Tables  3. 4-3. 6. 

Table  2.3  Ply  elastic  properties  for  IM7-8552. 


Property 

Value 

El  (GPa) 

171.42 

E2  (GPa) 

9.08 

Gi2  (GPa) 

5.29 

Vl2 

0.32 

Table  2.4  Ply  strengths  for  IM7-8552. 


Property 

Value  (MPa) 

Xt 

2326.2 

Xc 

1200.1 

Vud 

T 

62.3 

>C 

199.8 

^ud 

92.3 

Table  2.5  Calculated  in-situ  strengths  for  IM7-8552  (MPa). 


Ply  conhguration 

Vt 

Sl 

Thin  embedded  ply 

160.2 

130.2 

Thin  outer  ply 

101.4 

107.0 

Figures  2-2  and  2-3  show  respectively  the  held  variables  UVARMCl)  and  UVARM(3) 
of  the  0°  ply  for  an  applied  end  displacement  of  0.5mm. 
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2007 


Figure  2-2 


Field  variable  1  in  a  0° 


ply. 


Figure  2-3 


Field  variable  3  in  a  0° 


ply- 
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3.  UMAT  SUBROUTINE 


3.1  Overview 


The  second  part  of  this  report  is  a  reference  manual  for  the  use  of  an  ABAQUS® 
user  subroutine  UMAT  [1]  with  the  implementation  of  a  continuum  damage  model  for 
laminated  composites. 

The  plane-stress  continuum  damage  model  implemented  in  ABAQUS®  is  de¬ 
scribed  in  detail  in  references  [8]-[ll].  The  continuum  damage  model  is  dehned  in 
the  framework  of  the  thermodynamics  of  irreversible  processes.  Generally  speaking, 
the  formulation  of  the  continuum  damage  models  starts  by  the  dehnition  of  a  poten¬ 
tial  (for  example,  the  complementary  free  energy  per  unit  volume)  as  a  function  of  the 
damage  variables.  The  potential  is  the  basis  for  establishing  the  relation  between  the 
stress  and  the  strain  tensors.  For  the  complete  dehnition  of  the  constitutive  models 
it  is  also  necessary  to  dehne  the  damage  activation  functions,  i.e.  the  conditions  that 
lead  to  the  onset  of  inelastic  response,  and  the  damage  evolution  functions. 

The  present  continuum  damage  for  ABAQUS®  predicts  the  onset  and  accumu¬ 
lation  of  intralaminar  damage  mechanisms  (matrix  cracking  and  hber  fracture)  in 
laminated  composites  as  well  as  dual  structural  collapse  by  the  propagation  of  a 
macro-crack.  The  macro-crack  is  represented  by  a  line  of  localized  shell  or  continuum 
elements,  i.e.,  elements  where  the  constitutive  tangent  tensor  is  not  positive  dehnite. 

This  report  presents  the  basic  equations  required  for  the  dehnition  of  the  material 
properties  and  it  explains  how  to  dehne  the  model  in  ABAQUS®  standard.  An 
example  of  the  simulation  of  fracture  of  an  open-hole  carbon-epoxy  quasi-isotropic 
laminate  loaded  in  tension  is  also  described. 

The  full  details  of  the  development  and  validation  of  the  model  are  presented  in 
the  papers  shown  in  Appendix  B  and  Appendix  C: 


•  P.P.  Camanho,  P.  Maimi,  C.G.  Davila,  Prediction  of  size  effects  in  notched  lam¬ 
inates  using  continuum  damage  mechanics,  Composites  Science  and  Tech¬ 
nology,  67,  2715-2727,  2007.  This  paper  describes  the  application  of  the  model 
in  the  prediction  of  size  ehects  in  notched  composites  using  shell  elements. 

•  H.  Koerber,  P.P.  Gamanho,  Simulation  of  progressive  damage  in  bolted  compos¬ 
ite  joints,  Proceedings  of  the  13th  Enropean  Conference  on  Composite 
Materials  (ECCM-13),  Stockholm,  Sweden,  2008.  This  paper  presents  the 
application  of  the  model  in  the  prediction  of  the  mechanical  response  of  GFRP 
bolted  joints  using  solid  elements.  In  addition,  the  paper  compares  the  model 
proposed  here  with  the  ABAQUS®  built-in  damage  model  [1]. 
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3.2  Input  into  Abaqus  standard 


3.2.1  Shell  elements 


Material  properties 

The  material  properties  must  be  defined  by  the  user  in  the  j  obncone  .  inp  file  according 
to  the  following  example. 

** 

**  MATERIAL  #1 :  thin  embedded  ply 
** 

**  El,  E2,  G12,  vl2,  alphal , alpha2 ,  XT,  XPO 

**  XC,  YT,  YC,  ALPHAO,  SL  ,  SLud  ,  GIC_F1,  GIC_FE 

**  GIC_M,  GIIC_M,  GIC_FC,  GIIC_M~,  betal,  beta2  ,  DM,  Eta_viscous 

** 

^MATERIAL , NAME=IM7-8552-thin 

*USER  MATERIAL,  C0NSTANTS=24,  UNSYMM 

** 

171420.,  9080.  ,  5290.,  0.32  ,  -5.5E-6,  25.8E-6,  2323.5,  232.3, 

1200.1,  160.2  ,  199.8  ,  0.925,  130.2,  92.3,  31.5,  50.0, 

0.2774,  0.7879,  106.3  ,  1.3092,  0.000,  0.005,  0.0,  0.000 

^DENSITY 

1590E-6 

*DEPVAR 

15 


The  material  properties  defined  after  the  *USER  MATERIAL  command  are  shown  in 
Table  3.1. 


Table  3.1  Material  properties  in  the  jobname.inp  file. 


Line 

Column 

1 

2 

3 

4 

5 

6 

7 

8 

1 

El 

E2 

Gi2 

V12 

an 

0^22 

Xt 

Xpo 

2 

Xc 

Yt 

ttO 

Sl 

Sf 

G\p 

/^E 

^1+ 

3 

G2+ 

Gq 

Gi_ 

G2- 

/dll 

P22 

AM 

^visc 

Most  of  the  properties  required  are  standard  ply  properties.  However,  non¬ 
standard  material  properties  are  also  required.  Xpo  is  the  transition  strength  used 
for  the  definition  of  the  damage  evolution  law  for  longitudinal  tensile  failure. 
and  are  respectively  the  fracture  energies  per  unit  area  related  with  the  linear 
and  exponential  softening  laws  shown  on  Figure  3-1. 
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Figure  3-1  Damage  evolution  in  longitudinal  tension. 

The  parameter  ao  is  the  fracture  angle  of  a  ply  under  pure  transverse  compression. 
For  graphite/epoxy  and  glass/epoxy  materials  the  fracture  angle  ao  is  typically  equal 
to  53°.  is  the  parameter  used  in  the  viscous  regularization  procedure  used  to 
mitigate  convergence  difficulties. 

Initial  conditions 

The  model  requires  the  definition  of  the  initial  values  of  the  state  variables.  Therefore, 
the  j  obname .  inp  file  must  include  the  following  command: 

** 

**  Initialization  of  the  state  variables: 

** 

**  ELSET,  rfT,  rmT,  rfC,  rmC  ,  AfC,  AmT,  AmC,  dl,  d2,  d6, 

**  dr/dr,  gf,  gm,  FIfX,  ALEA 
** 

^INITIAL  CONDITIONS,  TYPE=S0LUTI0N 
OUT.PLT,  1.0,  1.0,  1.0,  1.0,  -1.0,  -1.0,  -1.0, 

0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  1.0 

where  0UT_PLT  represents  the  group  of  all  elements  whose  constitutive  model  is  de¬ 
fined  by  the  UMAT  subroutine.  The  state  variables  are  shown  in  Table  3.2. 
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Table  3.2  State  variables. 


STATEV(l) 

Damage  threshold  function 

STATEV(2) 

r2+ 

Damage  threshold  function 

STATEV(3) 

ri_ 

Damage  threshold  function 

STATEV(4) 

r2- 

Damage  threshold  function 

STATEV(5) 

Ai_ 

Adjustment  constant 

STATEV(6) 

A2+ 

Adjustment  constant 

STATEV(7) 

A2- 

Adjustment  constant 

STATEV(8) 

di 

Damage  variable 

STATEVO) 

d2 

Damage  variable 

STATEV(IO) 

de 

Damage  variable 

STATEV(ll) 

Derivative  required  for  the  viscous  regularization 

STATEV(12) 

9i 

Energy  dissipated 

STATEV(13) 

92  +  96 

Energy  dissipated 

STATEV(14) 

Damage  activation  functions 

STATEV(15) 

Ran 

Definition  of  random  properties 

The  state  variables  STATEV(l)  to  STATEV(14)  are  defined  in  detail  in  references 
[8]-[ll].  STATEV(15)  should  be  equal  to  one  if  a  random  field  of  material  properties  is 
not  required.  Setting  STATEV(15)=0  creates  a  random  field  for  the  ply  longitudinal 
strengths.  All  the  variables  can  be  post-processed  using  ABAQUS®  viewer: 

^ELEMENT  OUTPUT,  ELSET=DAMAGE_ELEMS 

1,  2,  3,  4 

SDV 


3.2.2  Continuum  elements 
Material  properties 

The  material  properties  must  be  defined  by  the  user  in  the  j  obname  .  inp  file  according 
to  the  following  example. 


26 


**  MATERIAL  #1 :  thin  embedded  ply 


** 

El  , 

E2 

G12 

vl2  ,  alphal. 

alpha2 , 

XT 

XPO, 

** 

XC 

YT 

YC 

ALPHAO  ,  SL 

SLud  , 

GIC.Fl, 

GIC.FE 

** 

GIC.M, 

GIIC.M, 

GIC_FC, 

GIIC.M-,  betal  , 

beta2  , 

DM 

Eta_viscous 

**  v23  ,  thickness 

** 

^MATERIAL , NAME=IM7-8552-thin 

*USER  MATERIAL,  C0NSTANTS=26 ,  UNSYMM 

** 

171420.,  9080.  ,  5290.,  0.32  ,  -5.5E-6,  25.8E-6,  2323.5,  232.3, 
1200.1,  160.2  ,  199.8  ,  0.925,  130.2,  92.3,  31.5,  50.0, 

0.2774,  0.7879,  106.3  ,  1.3092,  0.000,  0.005,  0.0,  0.000, 

0.52  ,0.125 

** 

^DENSITY 

1590E-6 

*DEPVAR 

15 


The  material  properties  defined  after  the  *USER  MATERIAL  command  are  shown  in 
Table  3.3. 


Table  3.3  Material  properties  in  the  jobname.inp  file. 


Line 

Column 

1 

2 

3 

4 

5 

6 

7 

8 

1 

El 

E2 

Gi2 

Vl2 

Oil 

0^22 

Xt 

Xpo 

2 

Xc 

Yt 

Oo 

Sl 

Sf 

Gt 

/^E 

^1+ 

3 

G2+ 

Gq 

Gi- 

G2- 

Pii 

P22 

AM 

^visc 

1^23 

t 

The  parameter  t  is  the  ply  thickness. 

Initial  conditions 

The  model  requires  the  definition  of  the  initial  values  of  the  state  variables.  Therefore, 
the  j  obname .  inp  file  must  include  the  following  command: 

** 

**  Initialization  of  the  state  variables: 

** 

**  ELSET,  rfT,  rmT,  rfC,  rmC  ,  AfC,  AmT,  AmC,  dl,  d2,  d6, 

**  dr/dr,  gf,  gm,  FIfX,  ALEA 
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^INITIAL  CONDITIONS,  TYPE=SOLUTION 
OUT.PLT,  1.0,  1.0,  1.0,  1.0,  -1.0,  -1.0,  -1.0, 
0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  1.0 


The  state  variables  are  shown  in  Table  3.2. 


3.2.3  Effect  of  Element  Size  and  Toughness  on  Degradation  Rates 

To  predict  damage  propagation,  it  is  necessary  to  establish  rates  of  material  degra¬ 
dation  that  are  energetically  consistent.  The  present  damage  model  relies  on  damage 
evolution  laws  that  take  into  account  the  toughness  of  the  material  in  each  damage 
mode  as  well  as  the  volume  of  material  represented  by  each  element  integration  point 
in  the  model  [8]-[ll].  However,  correct  degradation  rates  can  only  be  achieved  when 
the  hnite  element  mesh  is  sufficiently  rehned. 

When  an  element  is  too  large  for  a  given  material,  it  is  not  possible  to  achieve 
the  proper  degradation  without  a  snap-back  of  the  constitutive  model  [8]- [11].  Un¬ 
der  these  circumstances,  the  model  is  designed  to  reduce  the  material  strengths  as 
necessary  to  achieve  the  correct  energy  release  rates.  This  approach  allows  the  use  of 
coarser  meshes  for  damage  propagation.  However,  the  mesh  in  the  regions  of  damage 
initiation  should  be  sufficiently  rehned.  In  addition,  care  should  be  applied  to  ensure 
that  the  strengths  of  elements  in  coarse  regions  away  from  the  damage  areas  are  not 
reduced  so  much  as  to  cause  unrealistic  failures.  References  [8]- [11]  provide  techniques 
to  estimate  the  maximum  size  of  elements  in  the  regions  of  damage  initiation;  the 
maximum  size  for  the  hnite  element  for  each  damage  law  M  is: 

I*  <  M  =  1±,  2±,  6  (3.1) 

where  Em,  Gm  and  Xm  are  the  Young  modulus,  fracture  energies  and  strengths, 
respectively.  I*  is  the  characteristic  size  of  the  hnite  element  (ABAQUS®  CELENT 
parameter). 

When  the  strength  of  an  element  is  reduced  the  subroutine  writes  a  warning 
message  to  the  modelname.dat  hie  according  to  the  following  format: 

STRENGTH  REDUCTION  YT  =  101.229471616873  N.  ELEMENT  =  203 

where  YT  is  the  adjusted  value  of  the  transverse  tensile  strength. 

In  addition,  the  coarse  elements  whose  longitudinal  tensile  strength  and  transverse 
tensile  and  compressive  strengths  were  reduced  can  be  visualized  using  ABAQUS® 
viewer.  The  elements  whose  longitudinal  compressive  strength  {Xc)  was  reduced  have 
the  state  variable  STATEV(5)  equal  to  100.  The  elements  whose  transverse  tensile  and 
compressive  strengths,  Yr  and  Yc  respectively,  were  reduced  have  the  state  variable 
STATEV(6)  (Yr)  and  STATEV(7)  (Yc)  equal  to  100. 


3.2.4  Thermal  stresses 


The  constitutive  model  calculates  the  residual  thermal  stresses  that  result  from  the 
different  coefficients  of  thermal  expansion  in  the  longitudinal  and  transverse  direc¬ 
tions.  To  enable  the  calculation  of  the  thermal  stresses,  the  user  should  dehne  in  the 
jobname .  inp  hie  the  amplitude  of  the  thermal  step  as  follows: 

^AMPLITUDE,  NAME=AMPL,  DEFINITION=TABULAR 

The  residual  thermal  stresses  should  be  calculated  in  the  initial  step,  as  shown  in 
the  following  example: 

**  Step  1 :  thermal  step 

*STEP,  INC=10000,  UNSYMM=YES 
^STATIC 

0.5,  1.,  lD-7,  1. 

^TEMPERATURE , AMPLITUDE=AMPL 
GLOBAL,  -155. 

where  GLOBAL  represent  the  group  of  nodes  that  belong  to  the  elements  whose  consti¬ 
tutive  model  is  dehned  by  the  UMAT  subroutine,  and  -155  is  the  difference  between 
the  working  and  reference  temperatures. 

The  following  (mechanical)  steps  must  include  the  following  command: 

^TEMPERATURE 
GLOBAL,  -155. 


3.3  Example 


An  ABAQUS®  model  with  an  example  of  the  use  of  the  UMAT  subroutine  in  the 
strength  prediction  of  a  quasi-isotropic  Hexcel  IM7-8552  [90/0/ ±45] 3^  CFRP  laminate 
with  a  central  hole  and  loaded  in  tension  can  be  downloaded  from  the  following  URL: 

www.fe.up.pt/~pcamanho/oht3v2 . inp 

WWW . f e . up . pt/~pcamanho/model_open_hole_3_dy . inp 

www.fe.up.pt/~pcamanho/IM7-8552 . inp 

The  model  is  composed  of  ABAQUS®  S4  shell  elements  and  it  represents  a  speci¬ 
men  that  is  3mm  thick,  150mm  long,  36mm  wide  having  a  central  hole  with  a  diame¬ 
ter  of  6mm.  The  difference  between  the  working  and  stress-free  temperatures  used  to 
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calculate  the  residual  thermal  stresses  is  — 155°C.  An  implicit  dynamic  analysis  was 
performed  considering  a  loading  rate  of  2mm/min.  The  use  of  an  implicit  dynamic 
finite  element  model  enables  the  prediction  of  the  load  drop  that  occurs  when  the 
specimens  fail  catastrophically.  An  implicit  static  analysis  normally  fails  to  converge 
at  the  maximum  load. 

The  material  properties  used  are  shown  in  Tables  3. 4-3. 6. 

Table  3.4  Ply  elastic  properties  for  IM7-8552. 


Property 

Value 

El  (GPa) 

171.42 

E2  (GPa) 

9.08 

Gi2  (GPa) 

5.29 

V12 

0.32 

Table  3.5  Ply  strengths  for  IM7-8552. 


Property 

Value  (MPa) 

Xt 

2326.2 

Xc 

1200.1 

Vud 

It 

62.3 

199.8 

^ud 

92.3 

Table  3.6  Calculated  in-situ  strengths  for  IM7-8552  (MPa). 


Ply  conhguration 

Yt 

Sl 

Thin  embedded  ply 

160.2 

130.2 

Thin  outer  ply 

101.4 

107.0 

Table  3.7  Fracture  energies  for  transverse  fracture  for  IM7-8552  (kJ/m^). 


Property 

Value 

GI2+ 

0.2774 

Ge 

0.7879 

The  fracture  energy  G2-  is  calculated  as  G2-  =  Gq/  coscto  with  ao  =  53°. 

As  explained  in  section  3.2.3,  the  UMAT  subroutine  verihes  the  size  of  the  elements 
in  the  beginning  of  the  analysis:  if  the  size  of  one  element  is  large  enough  to  cause  a 
snap-back  of  the  constitutive  model,  the  strength  of  that  element  is  reduced.  There¬ 
fore,  the  mesh  should  be  rehned  in  the  locations  where  damage  initiation  is  likely 
to  take  place.  It  should  be  noted  that  a  strength  reduction  may  cause  difficulties  in 


30 


Table  3.8  Fracture  energies  for  longitudinal  fracture  for  IM7-8552  (kJ/m^). 


Property 

Value 

Gi+ 

81.5 

Gi_ 

106.3 

coarse  regions  of  the  model  away  from  the  damage  zones  where  the  strength  could  be 
reduced  enough  to  cause  premature  damage. 

To  overcome  this  difficulty,  it  is  suggested  to  assign  different  material  properties 
to  the  coarse  elements  located  in  regions  where  no  damage  takes  place.  The  strategy 
proposed  consists  in  increasing  the  fracture  toughness  of  these  elements  to  avoid  the 
strength  reduction.  For  example: 

****  MATERIAL  #4:  linear  elastic  material 
**** 

**  El,  E2,  G12,  vl2,  alphal,  alpha2,  XT  ,  XPO, 

**  XC,  YT,  YC,  ALPHAO  ,  SL,  SLud  ,  GIC.Fl,  GIC.FE 

**  GIC_M,  GIIC_M,  GIC_FC,  GIIC_M-,  betal  ,  beta2  ,  DM  ,  Eta_viscous 

** 

*MATERIAL,NAME=LE 

*USER  MATERIAL,  C0NSTANTS=24,  UNSYMM 
171420.,  9080.  ,  5290.,  0.32  ,  -5.5E-6,  25.8E-6,  2323.5,  232.3, 

1200.1,  113.3  ,  199.8  ,  0.925,  106.9,  92.3  ,  31.5e5,  50.0e5, 

0.2774e5,  0.7879e5,  106. 3e5  ,  1.3092e5,  0.000,  0.005,  0.0,  0.000 

** 

^DENSITY 

1590E-6 

*DEPVAR 

15 


Figure  3-2  highlights  in  red  the  region  modeled  with  coarse  elements  with  increased 
toughness,  which  implies  a  linear-elastic  response  for  these  elements.  Figure  3-3  shows 
the  region  modeled  with  the  actual  material  properties. 
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SPECIMEN  0HT3V2 
ODB:  oht3v2.odb 


Step:  Step-2 
Increment  400: 


ABAQUS/STANDARD  Version  6.5- 


Step  Time  =  33.37 


Dct  12  04:28:54  GMT  Daylight  Time  2006 


Figure  3-2  Coarse  mesh 


elements  with  toughness  increased. 


SPECIMEN  0HT3V2 
ODB:  oht3v2.odb  ABAQUS/STANDARD  Version  6.5- 

Jtep:  Step-2 

Increment  400:  Step  Time  =  33.37 


■ 
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Figure  3-3  Refined  mesh  -  real  material  properties. 


The  damage  variable  dg  in  the  0°  ply  at  the  onset  of  damage  localization  is  shown 
in  Figure  3-4. 
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Figure  3-4  Onset  of  damage  localization. 


The  predicted  load-displacement  relation  is  shown  in  Figure  3-5,  where  it  can  be 
observed  that  an  implicit  dynamic  analysis  predicts  the  load  drop  that  occurs  when 
the  specimen  fails  catastrophically. 
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4.  VUMAT  SUBROUTINE 


4.1  Overview 


The  third  part  of  this  report  is  a  reference  manual  for  the  use  of  an  ABAQUS®  user 
subroutine  VUMAT  [1]  with  the  implementation  of  the  continuum  damage  model  for 
laminated  composites  described  in  section  3  and  in  the  Appendix  B.  There  are  several 
relevant  loading  scenarios  for  which  an  explicit  Unite  element  code  is  more  appropriate. 
For  example,  in  dynamically  loaded  composite  structures  or  in  problems  with  multiple 
contact  surfaces.  In  addition,  explicit  formulations  can  provide  solutions  for  problems 
that  suffer  from  severe  convergence  difficulties  when  implicit  finite  element  codes  are 
used. 

Following  the  strategy  used  in  the  development  of  the  UMAT  subroutine,  the  code 
described  in  this  section  is  able  to  simulate  the  mechanical  response  of  composite 
structures  using  both  shell  and  continuum  elements. 

This  section  includes  an  example  of  the  simulation  of  a  low-velocity  impact  in 
a  composite  laminate,  using  the  VUMAT  subroutine  developed  in  this  project,  an  ad¬ 
ditional  VUMAT  subroutine  that  simulates  delamination  onset  and  propagation  (the 
development  of  this  second  VUMAT  subroutine  was  planned  for  the  second  year  of  this 
project). 


4.2  Input  into  Abaqus  explicit 

4.2.1  Shell  elements 
Material  properties 

The  material  properties  must  be  defined  by  the  user  in  the  j  obname  .  inp  file  according 
to  the  following  example. 

** 

**  USER  PLY  MATERIAL 
** 

**  El  ,  E2  ,  G12  ,  vl2  ,  alphal,  alpha2,  XT  ,  XPO, 

**  XC  ,  YT  ,  YC  ,  ALPHAO  ,  SL  ,  SLud  ,  GIC.Fl,  GIC.FE 
**  GIC_M,  GIIC_M,  GIC_FC,  GIIC_M-,  betal  ,  beta2  ,  DM  ,  Eta_viscous 
** 

^MATERIAL ,  NAME=IM7-8552-Dainage 
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*USER  MATERIAL,  C0NSTANTS=25 

171420.,  9080.,  5290.,  0.32  ,  -5.5E-6,  25.8E-6,  2323.5,  232.3, 

1200.1,  132.7,  199.8,  0.925,  117.1,  92.3,  81.5,  10.0, 

1.0,  2.0,  106.3,  1.31,  0.000,  0.005,  0.0,  0.0, 

^Density 
1.59e-9 

*DEPVAR,DELETE=17 

17 

1,  rfT,  "Fiber  tension  internal  variable" 

2,  rmT,  "Matrix  tension  internal  variable" 

3,  rfC,  "Fiber  compr  internal  variable" 

4,  rmC,  "Matrix  compr  internal  variable" 

5,  AfC,  "Fiber  tension  adjustment  parameter" 

6,  AmT,  "Matrix  tension  adjustment  parameter" 

7,  AmC,  "Matrix  compr  adjustment  parameter" 

8,  dl,  "Damage  variable,  direction  11" 

9,  d2,  "Damage  variable,  direction  22" 

10,  d6,  "Damage  variable,  direction  12" 

11,  gf,  "Fiber  dissipated  energy" 

12,  gm,  "Matrix  dissipated  energy" 

13,  Ell,  "Direct  strain,  direction  11" 

14,  E22,  "Direct  strain,  direction  22" 

15,  E33,  "Direct  strain,  direction  33" 

16,  E12,  "Shear  strain,  direction  12" 

17,  STATUS,  "Status  of  the  element" 

The  STATUS  variable  defines  the  status  of  an  element:  if  STATUS=1  the  element 
is  active,  and  if  STATUS=0  the  element  has  been  deleted.  The  criterion  implemented 
to  delete  an  element  from  the  mesh  is  based  on  the  value  of  the  damage  variable 
associated  with  failure  mechanisms  in  the  longitudinal  direction;  if  di  >  0.999  the 
element  is  deleted  from  the  mesh.  This  procedure  prevents  the  severe  reduction  on  the 
stable  time  increment  that  results  from  highly  distorted,  damaged,  finite  elements. 

The  material  properties  defined  after  the  *USER  MATERIAL  command  are  shown 
in  Table  3.1. 

Initial  conditions 

The  model  requires  the  definition  of  the  initial  values  of  the  state  variables.  Therefore, 
the  j  obname .  inp  file  must  include  the  following  command: 
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**  shell  elements 

**  ELSET,  rfT,  rmT,  rfC,  rmC,  AfC,  AmT,  AmC, 

**  dl,  d2,  d6,  gf,  gm,  ell,  e22,  e33, 

**  el2,  STATUS 
** 

^Initial  Conditions,  Type=Solution 

<elset>,  1.0,  1.0,  1.0,  1.0,  -1.0,  -1.0,  -1.0, 

0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0, 

0.0,  1.0 

where  <elset>  represents  the  group  of  all  elements  whose  constitutive  model  is  de¬ 
fined  by  the  VUMAT  subroutine. 

The  state  variables  used  are  shown  in  Table  4.1. 


Table  A 

kl  State  variables. 

STATEV(l) 

Damage  threshold  function 

STATEV(2) 

Damage  threshold  function 

STATEV(3) 

ri_ 

Damage  threshold  function 

STATEV(4) 

r2- 

Damage  threshold  function 

STATEV(5) 

Ai_ 

Adjustment  constant 

STATEV(6) 

^2+ 

Adjustment  constant 

STATEV(7) 

^2- 

Adjustment  constant 

STATEV(8) 

dl 

Damage  variable 

STATEVO) 

d2 

Damage  variable 

STATEV(IO) 

de 

Damage  variable 

STATEV(ll) 

9i 

Energy  dissipated 

STATEV(12) 

92  +  96 

Energy  dissipated 

STATEV(13) 

^11 

11-component  of  the  strain  tensor 

STATEV(14) 

^22 

22-component  of  the  strain  tensor 

STATEV(15) 

^33 

33-component  of  the  strain  tensor 

STATEVde) 

^12 

12-component  of  the  strain  tensor 

STATEV(17) 

STATUS 

Status  of  the  element 

4.2.2  Continuum  elements 


Material  properties 

The  material  properties  must  be  defined  by  the  user  in  the  j  obname  .  inp  hie  according 
to  the  following  example. 
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**  MATERIAL  #1 :  thin  embedded  ply 
** 

**  El  ,  E2  ,  G12  ,  vl2  ,  alphal,  alpha2,  XT  ,  XPO, 

**  XC  ,  YT  ,  YC  ,  ALPHAO  ,  SL  ,  SLud  ,  GIC.Fl,  GIC.FE 
**  GIC_M,  GIIC_M,  GIC_FC,  GIIC_M-,  betal  ,  beta2  ,  DM  ,  Eta_viscous 
**  v23  ,  thickness 

** 

^MATERIAL , NAME=IM7-8552-thin 

*USER  MATERIAL,  C0NSTANTS=26 ,  UNSYMM 

** 

171420.,  9080.  ,  5290.,  0.32  ,  -5.5E-6,  25.8E-6,  2323.5,  232.3, 

1200.1,  160.2  ,  199.8  ,  0.925,  130.2,  92.3,  31.5,  50.0, 

0.2774,  0.7879,  106.3  ,  1.3092,  0.000,  0.005,  0.0,  0.000, 

0.52  ,0.125 

^DENSITY 

1590E-6 

*DEPVAR,DELETE=20 

20 

1,  rfT,  "Fiber  tension  internal  variable" 

2,  rmT,  "Matrix  tension  internal  variable" 

3,  rfC,  "Fiber  compr  internal  variable" 

4,  rmC,  "Matrix  compr  internal  variable" 

5,  AfC,  "Fiber  tension  adjustment  parameter" 

6,  AmT,  "Matrix  tension  adjustment  parameter" 

7,  AmC,  "Matrix  compr  adjustment  parameter" 

8,  dl,  "Damage  variable,  direction  11" 

9,  d2,  "Damage  variable,  direction  22" 

10,  d6,  "Damage  variable,  direction  12" 

11,  gf,  "Fiber  dissipated  energy" 

12,  gm,  "Matrix  dissipated  energy" 

13,  Ell,  "Direct  strain,  direction  11" 

14,  E22,  "Direct  strain,  direction  22" 

15,  E33,  "Direct  strain,  direction  33" 

16,  E12,  "Shear  strain,  direction  12" 

17,  E13,  "Shear  strain,  direction  13" 

18,  E23,  "Shear  strain,  direction  23" 

19,  d3,  "Damage  variable  direction  33" 

20,  STATUS,  "Status  of  the  element" 


The  material  properties  defined  after  the  *USER  MATERIAL  command  are  shown  in 
Table  3.3. 

The  state  variables  used  are  shown  in  Table  4.2. 
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Table  A 

l■.2  State  variables. 

STATEV(l) 

Damage  threshold  function 

STATEV(2) 

Damage  threshold  function 

STATEV(3) 

ri_ 

Damage  threshold  function 

STATEV(4) 

r2- 

Damage  threshold  function 

STATEV(5) 

^1- 

Adjustment  constant 

STATEV(6) 

^2+ 

Adjustment  constant 

STATEV(7) 

^2- 

Adjustment  constant 

STATEV(8) 

dl 

Damage  variable 

STATEVO) 

d2 

Damage  variable 

STATEV(IO) 

de 

Damage  variable 

STATEV(ll) 

9i 

Energy  dissipated 

STATEV(12) 

92  +  96 

Energy  dissipated 

STATEV(13) 

^11 

11-component  of  the  strain  tensor 

STATEV(14) 

^22 

22-component  of  the  strain  tensor 

STATEV(15) 

^33 

33-component  of  the  strain  tensor 

STATEVde) 

^12 

12-component  of  the  strain  tensor 

STATEV(17) 

£^13 

13-component  of  the  strain  tensor 

STATEVdS) 

^23 

23-component  of  the  strain  tensor 

STATEVd9) 

d3 

Damage  variable 

STATEV(20) 

STATUS 

Status  of  the  element 

Initial  conditions 

The  model  requires  the  dehnition  of  the  initial  values  of  the  state  variables.  Therefore, 
the  j  obname .  inp  hie  must  include  the  following  command; 

**  3D  elements 

**  ELSET,  rfT,  rmT,  rfC,  rmC,  AfC,  AmT,  AmC, 

**  dl,  d2,  d6,  gf,  gm,  ell,  e22,  e33, 

**  el2,  el3,  e23,  STATUS 
** 

^Initial  Conditions,  Type=Solution 

<elset>,  1.0,  1.0,  1.0,  1.0,  -1.0,  -1.0,  -1.0, 

0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0, 

0.0,  0.0,  0.0,  0.0,  1.0 
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4.3  Example 


An  ABAQUS®  model  with  an  example  of  the  use  of  the  VUMAT  subroutine  in  the 
strength  prediction  of  a  Hexcel  IM7-8552  [idSs]  CFRP  laminate  under  a  low- velocity 
impact  can  be  downloaded  from  the  following  URL: 

www.fe.up.pt/~pcamanho/Impact_Test . inp 

WWW . f e . up . pt/~pcamanho/Laminate_Solid_Test . id 

The  model  of  the  composite  plies  is  created  using  ABAQUS®  C3D8  solid  elements 
and  it  represents  a  2mm  thick,  lOOmmx  100mm  square  specimen,  as  shown  in  Figure 
4-1. 


Figure  4-1  Specimen  and  impactor. 

The  impact  results  from  the  contact  of  the  specimen  with  a  semi-spherical  rigid 
body  with  a  diameter  of  16mm,  mass  of  1kg,  initial  velocity  of  4m/s,  corresponding 
to  an  impact  energy  of  8J.  The  composite  specimen  is  clamped  along  all  edges.  The 
material  properties  used  are  shown  in  Tables  3. 4-3. 6.  In  addition  to  the  the  simulation 
of  ply  failure  mechanisms  by  means  of  the  subroutine  VUMAT,  the  separation  between 
the  -|-45  and  -45  plies  (delamination)  is  simulated  using  another  VUMAT  subroutine 
where  a  cohesive  formulation  previously  proposed  is  implemented  in  ABAQUS®  solid 
cohesive  elements  [12]. 
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Figure  4-2  shows  the  predicted  delaminated  region  at  the  interface  between  the 
-|-45  and  -45  plies.  The  damage  variable  dg  in  the  —45°,  corresponding  to  the  back-face 
of  the  laminate,  is  shown  in  Figure  4-3. 
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Figure  4-3  Predicted  transverse  matrix  cracking  in  the  laminate  back-face. 
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Figure  4-4  shows  the  damage  on  the  top  (+45)  ply,  at  the  contact  region  between 
the  impactor  and  the  laminate. 


Figure  4-4  Predicted  transverse  matrix  cracking  in  the  contact  region  between  the 
impactor  and  laminate. 


The  preliminary  results  obtained  in  this  example  indicate  that  the  model  devel¬ 
oped  is  able  to  simulate  the  interaction  between  the  failure  mechanisms.  In  addition, 
the  implementation  of  the  material  model  in  ABAQUS®  explicit  renders  the  solution 
of  complex  dynamic  problems  involving  contact  and  failure  possible. 
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5.  CONCLUSIONS 


A  new  ABAQUS®  UVARM  subroutine  with  the  computational  implementation  of  the 
LaRC  failure  criteria  was  developed  in  the  hrst  phase  of  this  project.  This  subroutine 
predicts  hrst  ply  failure,  and  it  may  be  used  for  the  preliminary  analysis  of  composite 
structures,  for  the  identihcation  of  critical  regions  of  such  structures,  and  in  the  design 
of  structures  where  no  type  of  damage  is  acceptable  (e.g.  cryogenic  fuel  tanks). 

A  continuum  damage  model  able  to  predict  the  onset  and  propagation  of  in¬ 
tralaminar  failure  mechanisms  was  developed  and  implemented  in  both  UMAT  and 
VUMAT  ABAQUS®  subroutines.  The  model  may  be  used  in  both  shell  and  three- 
dimensional  continuum  elements  using  ABAQUS®  implicit  (UMAT)  and  ABAQUS® 
explicit  (VUMAT).  The  computational  models  were  implemented  according  to  the  re¬ 
quirements  establishes  in  the  beginning  of  the  project: 

•  Accurate  prediction  of  damage  onset.  The  failure  criteria  implemented  is  able 
to  represent  the  following  characteristics  of  the  mechanical  behavior  of  lami¬ 
nated  composite  materials:  i)  in-situ  effects,  i.e.  the  effective  increase  of  the 
transverse  tensile  and  in-plane  shear  strengths  of  a  ply  when  it  is  embedded  in 
a  multidirectional  laminate;  ii)  the  beneficial  effect  of  transverse  compression 
on  the  apparent  shear  strength  of  a  ply;  and  hi)  the  effect  of  the  shear  stresses 
on  fiber  kinking. 

•  Crack  closure  under  load  reversal.  The  continuum  damage  model  implements 
an  unilateral  representation  of  cracks,  allowing  the  load  path  continuity  to  be 
recovered  when  cracks  close  under  compressive  loads. 

•  Residual  thermal  stresses.  The  constitutive  model  represents  the  effects  of  the 
residual  thermal  stresses  in  the  plies  of  a  multidirectional  laminate. 

•  Standard  material  properties.  The  majority  of  the  material  properties  required 
by  the  model  can  be  obtained  from  standard  test  methods. 

•  Ply-level  properties.  The  model  uses  ply  properties,  thus  avoiding  the  need  to 
test  laminates  every  time  the  lay-up  or  stacking  sequence  is  modified. 

•  Regularization  of  energy  dissipated.  The  model  avoids  mesh  dependency  prob¬ 
lems  and  assures  the  objectivity  of  the  numerical  solution  by  accounting  for  the 
toughnesses  of  the  material  in  each  damage  mode  as  well  as  the  energy  dissi¬ 
pated  by  damage  at  a  material  integration  point.  In  addition,  procedures  to 
rapidly  assess  adequacy  of  the  mesh  resolution  and  to  provide  corrective  means 
when  maximum  mesh  size  requirements  cannot  be  met  were  proposed. 

•  Explicitly  integrated  constitutive  model.  The  model  does  not  require  iterations 
to  solve  the  constitutive  equations,  being  therefore  suitable  to  be  used  in  large 
scale  computations. 
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•  Fast  convergence  rate  of  numerical  solution.  The  model  includes  stabilization 
methods  that  mitigate  the  problem  of  slow  convergence  rate  of  models  involving 
strain-softening. 

It  was  also  shown  that  the  continuum  damage  model  used  in  shell  elements  is 
able  to  predict  the  effect  of  size  on  the  strength  of  quasi-isotropic  CFRP  laminates. 
In  addition,  the  preliminary  validation  example  of  a  low-velocity  impact  load  on  a 
CFRP  laminate  indicates  that  the  combination  of  a  ply-based  damage  model  and  a 
cohesive  formulation  for  the  simulation  of  delamination  is  a  strategy  that  may  capture 
the  interaction  between  these  different  failure  mechanisms. 

Future  work  should  address  the  issues  related  to  the  mesh-dependent  directional¬ 
ity  of  crack  propagation  that  is  often  observed  in  continuum  damage  models.  This 
problem  may  be  mitigated  by  using  an  improved  kinematic  representation  of  the  fail¬ 
ure  mechanisms,  such  as  the  transversely  isotropic  damage  model  presented  in  [13] 
and  presented  in  Appendix  D  (the  full  development  of  this  model  was  planned  for  the 
second  year  of  the  project). 
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Appendix  A:  LaRC03-UVARM-vl  Fortran  subroutine 


SUBROUTINE  UVARM (UVAR .DIRECT , T , TIME , DTIME , CMNAME , ORNAME , 

1  NUVARM ,  NOEL ,  NPT ,  LAYER ,  KSPT ,  KSTEP ,  KING ,  NDI ,  NSHR ,  COORD , 

2  JMAC , JMATYP , MATLAYO , LACCFLA ) 

C 

INCLUDE  ' ABA_PARAM . INC ’ 

C 

common/crdf Ig/Irdf Ig 
C 

CHARACTER*80  CMNAME , ORNAME , CMNAMEl 
CHARACTER'S  FLGRAY(15) 

CHARACTER  xoutdir*255,  xfname*80 
CHARACTER  dmkname*255 ,  FNAMEX*80 
DIMENSION  UVAR(*) ,DIRECT(3,3) ,T(3,3) ,TIME(2) 

DIMENSION  ARRAY(15) , JARRAY(15) , JMAC(*) , JMATYPC*) ,C00RD(*) 

C 

DIMENSION  stress (6) 

C  The  dimensions  of  the  variables  FLGRAY,  ARRAY  and  JARRAY 

C  must  be  set  equal  to  or  greater  than  15. 

double  precision  aIphao,aIphamem(l) ,psimem(l) ,thetamem(l) , 

1  omega ( 1 ) , lambda ( 1 ) , 

2  fmat(l) ,fkink(l) ,fft(l) ,epsmato(l) , 

3  sigmato(l) ,epskinko(l) ,sigkinko(l) ,epsfto(l) , 

4  sigl(l) ,sig2(l) ,sig3(l) ,taul2(l) ,tau23(l) ,tau31(l) , 

5  epsl(l) ,eps2(l) ,eps3(l) ,epsl2(l) ,eps23(l) ,eps31(l) , 

6  sl2,s23,fio,beta 
C 

integer  1ft, lit 
C 

pi=dacos (-1 . dO) 
degtorad=pi/180 . dO 
C 

do  i=l,nuvarm 
uvar(i)  =  O.dO 
enddo 
C 


C  Open  and  read  input  file  with  material  properties: 

C  directory/jobname .mt 


Ixfname  =  0 
Ixoutdir  =  0 
xfname  ='  ’ 
xoutdir  =’  ’ 
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call  get jobnameCxf name, Ixf name)  !  input  file  name 

call  getoutdir(xoutdir,lxoutdir)  !  output  directory 

if (Irdf Ig.ne . 1)  then 

fnamex=dmkname (xfname (1 : Ixf name) ,xoutdir (1 : Ixoutdir) , ’ .mt ’ ) 
open (unit=17 , f ile=f namex , status= ’ old ’ ) 

Irdflg  =  1 
endif 

read  (17,*) 
read  (17,*)  klarc 

CMNAME1= ' **dummy_name** ’ 

do  while (CMNAMEl .NE.CMNAME)  !  search  for  material  type 

read  (17,*) 
read  (17,*)  CMNAMEl 
if (CMNAMEl. EQ.CMNAME)  then 
read  (17,*) 

read  (17,*)  yml,  ym2,  ym3,  nu21,  nu31,  nu32 
read  (17,*) 

read  (17,*)  gl2,  g23,  g31,  xt,  xc,  yt,  yc,  sl2 
read  (17,*) 

read  (17,*)  alphao,  beta,  g,  slis 
else 

do  i=l,6 
read(17, *) 
enddo 
endif 
enddo 

rewind  17 


C  Compute  derived  material  properties 


alphao=alphao*degtorad 

st=yc*dcos (alphao) * (dsin (alphao) +dcos (alphao) /dtan(2 . dO*alphao) ) 

s23=st 

sl=sl2 

etat=-l . d0/dtan(2 .dO*alphao) 

etal=-sl2*dcos (2*alphao) / (yc*dcos (alphao) *dcos (alphao) ) 


C  Read  stress  tensor  from  current  increment 


CALL  GETVRM (’S’, ARRAY , JARRAY , FLGRAY , JRCD , JMAC , JMATYP , MATLAYO , 
1  LACCFLA) 
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!  LaRC03  failure  criteria 


C 


C 


if (klarc . eq. 
stress(l)  = 
stress (2)  = 
stress (3)  = 
stress (4)  = 
stress (5)  = 
stress (6)  = 


3)  then 

■  array (1) 

■  array (2) 

■  array (3) 

■  array (4) 

■  array (5) 

■  array (6) 


call  larc03(stress(l) ,stress(2) ,stress(3) ,stress(4) , 

1  stresses) , stresses) ,XT,XC,YT, YC, 

2  SL , SLIS , ST , G , G12 , ETAL , ETAT , NDIM , UVAR , ANGLES , 

3  NOUT.NUVARM) 
endif 


* 

*  End  of  main  program 

* 


RETURN 

END 

*  <<<<<<<<<<<<<<<<<<<<<<<<  SUBROUTINE  LARGOS  >>»>»>>»»»>>»»»>>  * 

*  * 

*  LaRCOS  failure  criteria  * 

*  * 


*  <«<<«««<<«««<<«««<<«<«>>>»»»>>»>»>>»»»>>»»»>>  * 
SUBROUTINE  LaRCOS eSll , S22 , S33 , S12 , S13 , S23 , XT , XC , YT , YC , 

1  SL , SL_IS , ST , G , G12 , ETA_L , ETA_T , NDIM , FI , ANGLES , 

2  NOUT,NUVARM) 

C 

IMPLICIT  NONE 
C 

DOUBLE  PRECISION  Sll , S22 , S33 , S12 , S13 , S23 ,XT , XC , YT, YC , 

1  SL,SL_IS,ST,G,G12,ETA_L,ETA_T,Fie*) .ANGLES, 

2  PI , SI 1_M , S22_M , S33_M , S12_M , S13_M , S23_M , FLARC03 , 

3  FMCCAULEY, ALPHA, FIP ex) 

C 

INTEGER  NDIM,NOUT,I,NUVARM 
C 

PI  =  DACOSe-l.DO) 

C 

do  i=l,7 
fipei)  =  O.dO 
enddo 

C - 

* 

*  Transverse  ematrix) 

* 
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IF(S22.GT.0.D0)  THEN  !  matrix  tension 

IF(S11.LT.0.D0.AND.DABS(S11) .LT.XC/2)  THEN 

CALL  RDTATE_PHI (SL , XC , ETA_L , SI 1 , S22 , S12 , 0 . DO , 0 . DO , 

1  S11_M,S22_M,S12_M,S13_M,S23_M,G12,NDIM) 

FIP(l)  =  (1-G)*S22_M/YT+G*S22_M/YT*S22_M/YT+S12_M/SL_IS* 

1  S12_M/SL_IS 

ELSE 

FIP(2)  =  (1-G)*S22/YT+G*S22/YT*S22/YT+S12/SL_IS*S12/SL_IS 
ENDIF 
C 

ELSE  !  matrix  compression 

IF(Sll.GE.-YC)  THEN 

FIP(3)  =  FLaRC03(ALPHA,S22,S12,ETA_L,ETA_T,SL_IS,ST,PI) 
ELSE 

CALL  RDTATE_PHI (SL , XC , ETA_L , SI 1 , S22 , S12 , 0 . DO , 0 . DO , 

1  S11_M,S22_M,S12_M,S13_M,S23_M,G12,NDIM) 

FIP(4)  =  FLaRC03(ALPHA,S22_M,S12_M,ETA_L,ETA_T,SL_IS,ST,PI) 
ENDIF 
ENDIF 

C - 

* 

*  Longitudinal  (fibre) 

*  C 


IF(Sll.GE.O.DO)  THEN  !  fibre  tension 

FIP(5)  =  Sll/XT 
C 

ELSE  !  fibre  compression 

CALL  ROTATE_PHI (SL , XC , ETA_L , SI 1 , S22 , S12 , 0 . DO , 0 . DO , 

1  S11_M,S22_M,S12_M,S13_M,S23_M,G12,NDIM) 

IF(S22_M.LT.0.D0)  THEN 

FIP(6)  =  FMcCAULEY((DABS(S12_M)+ETA_L*S22_M)/SL_IS)  !  LaRC#4 
ELSEIF  (DABS(Sll) .GE.XC/2.D0)  THEN 
FIP(7)  =  (1-G)*S22_M/YT+G*S22_M/YT*S22_M/YT+S12_M/SL_IS 
1  *S12_M/SL_IS 

ENDIF 

ENDIF 

C 

FI(1)  =  MAX(FIP(1) ,FIP(2))  !  Transverse  with  S22>0 

C 

FI (2)  =  MAX(FIP(3) ,FIP(4))  !  Transverse  with  S22<0 

C 

FI (3)  =  FIP(5)  !  Longitudinal  with  S11>0 

C 

FI(4)  =  MAX(FIP(6)  ,FIP(7))  !  Longitudinal  with  SIKO 

C 
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* 

* 

* 

* 

* 

c 

c 


RETURN 

END 

<<<<<<<<<<<<<<<<<<<<<<<<<  FUNCTION  FLaRCOS  >>»>»>>»»»>>»»»>>  * 

* 

MATRIX  COMPRESSION  FAILURE  CRITERION  (LaRCOS)  * 

* 

<«<<«««<<«««<<«««<<«»»>>»»»>>»>»>>»»»>>»»»>>  * 
REALMS  FUNCTION  FLaRC03(ALPHA,S22,S12,ETAL,ETAT,SL_IS,ST,PI) 

IMPLICIT  NONE 

DOUBLE  PRECISION  ALPHA, S22 , S12 ,ETAL , ETAT, SL_IS , ST, PI , THETA , 

1  TAUT_EFF,  FMCCAULEY,TAUL_EFF,FIC,ALPHAl,tauI,taut 

INTEGER  I 

Cycle  over  possible  fracture  angles 
FLaRC03=0 . dO 

DO  i=0,56  !  Determination  of  the  fracture  angle 
ALPHAl  =  i*PI/180.D0 

IF(ALPHA1.EQ.0.D0.0R.S22.EQ.0.D0)  THEN  !  Avoids  divisions  by  zero 
THETA  =  PI/2. DO 
ELSE 

THETA  =  DATAN(-dabs(S12)/(S22*DSIN(ALPHAl))) 

ENDIF 

TAUT_EFF  =  FMcCAULEY(-S22*dcos(alphal)*(dsin(alphal)-etat* 

1  dcos (alphal) *dcos (theta) ) ) 

TAUL_EFF  =  FMcCAULEY(dcos(alphal)*(dabs(sl2)+etal*s22* 

1  dcos (alphal) *dsin(theta) ) ) 


FIC 


1 


(TAUT_EFF/ST) * (TAUT_EFF/ST) + 
(TAUL_EFF/SL_IS) * (TAUL_EFF/SL_IS) 


FLaRCOS  =  max(FLaRC03,FIC) 
ENDDO 


RETURN 

END 

*  <<<<<<<<<<<<<<<<<<<<<<  SUBROUTINE  ROTATE_PHI>»>»>>»»»>>»»»>>  * 

*  * 

*  ROTATION  OF  STRESSES  TO  THE  MISALIGNMENT  COORDINATE  FRAME  * 

*  * 


*  <«<<«««<<«««<<«««<<«<«>>>»»»>>»>»>>»»»>>»»»>>  * 
SUBROUTINE  ROTATE_PHI ( SL , XC , ETAL ,S11,S22,S12,S13,S23, 

1  S11T,S22T,S12T,S13T,S23T,G12,NDI) 

C 

IMPLICIT  NONE 
C 
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DOUBLE  PRECISION  SL , XC ,ETAL , Sll , S22 , S12 , S13 , S23 , 

1  SllT,S22T,S12T,S13T,S23T,G12,aa,cc,phiC, 

2  PHI,sqr,phiO,cp,ss,c2,s2 

INTEGER  NDI 

Calculate  fiber  misalignment  angle  (Linear  shear  law  and 
small  angle  approximations) 
cc  =  dABS(SL/XC) 
aa  =  cc+ETAL 

sqr  =  dsqrt (1 . dO-4 . OdO*aa*cc) 

phiC  =  datan( (1 .dO-sqr) / (2 . OdO*aa) )  !  select  smallest  root 

phiO  =  (dabs(S12)+(G12-XC)*phiC)/(G12+Sll-S22) 

cp  =  dcos(phiO) 

ss  =  dsin(phiO) 

c2  =  cp*cp 

s2  =  ss*ss 

C  Calculate  stresses  in  misalignment  coordinate  frame  C 
SllT  =  Sll*c2+S22*s2+2.0d0*cp*ss*DABS(S12) 

S22T  =  Sll*s2+S22*c2-2.0d0*cp*ss*DABS(S12) 

S12T  =  -ss*cp=t=Sll+ss*cp*S22+(c2-s2)*DABS(S12) 

RETURN 

END 

<<<<<<<<<<<<<<<<<<<<<<<  FUNCTION  FMcCAULEY  >>»>»>>»»»>>»»»>>  * 

* 

McCauley  operator  * 

* 

<«<<«««<<«««<<«««<<«»»>>»»»>>»>»>>»»»>>»»»>>  * 
REAL*8  FUNCTION  FMcCAULEY (X) 

IMPLICIT  NONE 

DOUBLE  PRECISION  X 

IF(X.LE.O.DO)  THEN 
FMcCAULEY  =  O.DO 
ELSE 

FMcCAULEY  =  X 
ENDIF 

RETURN 

END 

<<<<<<<<<<<<<<<<<<<<<<<<  FUNCTION  DMKNAME  >>»>»>>»»»>>»»»>>  * 


Compose  a  filename  directory/jobname.exten 
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*  <«<<«««««««<<«««««»»>>»»»>>»>»>>»»»>>»»»>>  * 

character* (*)  function  dmkname(fname,dname,exten) 

C 

character* (*)  f name, dname, ext en 

C  C  fname  I  jobname  C  dname  I  directory  C  exten  I 
extension  C  dmkname  0  directory/jobname . exten  C 
Itot  =  len(fname) 

If  =  0 

do  kl  =  Itot, 2,-1 

if  (If . eq. 0 . and. fname (kl :kl) .ne . ’  ’)  If  =  kl 
end  do 
C 

Itot  =  len (dname) 

Id  =  0 

do  kl  =  Itot, 2,-1 

if  (Id. eq. 0 . and. dname (kl :kl) .ne . ’  ’)  Id  =  kl 
end  do 
C 

Itot  =  len(exten) 
le  =  0 

do  kl  =  Itot, 2,-1 

if  (le . eq. 0 . and. exten(kl :kl) .ne . ’  ’)  le  =  kl 
end  do 
C 

if  ((If  +  Id  +  le)  .le.  len(dmkname) )  then 
dmkname  =  dname (1 : Id) // ’/’//fname (1 : If ) 

Itot  =  Id  +  If  +  1 
if  (  le.gt.O)  then 

dmkname  =  dmkname  (1 :  Itot) //extend  :  le) 
end  if 
end  if 
C 

return 

end 


C  ====  end  of  program  ==== 
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Appendix  B:  Paper  published  in  Composites  Science 
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Abstract 

This  paper  examines  the  use  of  a  continuum  damage  model  to  predict  strength  and  size  effects  in  notched  carhon-epoxy  laminates. 
The  effects  of  size  and  the  development  of  a  fracture  process  zone  before  final  failure  are  identified  in  an  experimental  program.  The 
continuum  damage  model  is  described  and  the  resulting  predictions  of  size  effects  are  compared  with  alternative  approaches:  the  point 
stress  and  the  inherent  flaw  models,  the  Linear  Elastic  Fracture  Mechanics  approach,  and  the  strength  of  materials  approach.  The  results 
indicate  that  the  continuum  damage  model  is  the  most  accurate  technique  to  predict  size  effects  in  composites.  Furthermore,  the  contin¬ 
uum  damage  model  does  not  require  any  calibration  and  it  is  applicable  to  general  geometries  and  boundary  conditions. 

©  2007  Elsevier  Ltd.  All  rights  reserved. 
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1.  Introduction 

The  introduction  of  advanced  composite  materials  in 
new  applications  relies  on  the  development  of  accurate 
analytical  and  computational  tools  that  are  able  to  predict 
the  thermo-mechanical  response  of  composites  under  gen¬ 
eral  loading  conditions  and  geometries.  In  the  absence  of 
accurate  analytical  models,  the  design  process  has  to  rely 
on  costly  matrices  of  mechanical  tests  based  on  large  num¬ 
bers  of  test  specimens  [1]  and  empirical  knockdown  factors 
[2]. 

The  prediction  of  ultimate  strength  remains  the  main 
challenge  in  the  simulation  of  the  mechanical  response  of 
composite  materials  [3].  The  simulation  of  size  effects  on 
the  strength  of  composites  is  of  particular  interest  and  rel¬ 
evance  [4—8]:  reliable  analytical  and  numerical  models  must 
represent  the  decrease  of  the  ultimate  strength  when  the 
structural  dimensions  increase  [9]. 
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Size  effects  in  laminated  composites  occur  at  different 
material  and  structural  levels.  At  the  meso-mechanical 
level,  it  is  observed  that  the  transverse  tensile  and  in-plane 
shear  strengths  of  a  ply  constrained  by  sublaminates 
depend  on  the  ply  thickness  [10].  This  size  effect  is  normally 
called  the  “in  situ”  effect  and  can  be  accounted  for  in  the 
prediction  of  matrix  cracking  onset  using  the  “in  situ” 
strengths  in  appropriate  failure  criteria.  The  “in  situ” 
strengths  can  be  calculated  from  analytical  closed-fomi 
solutions  using  ply  elastic  properties  and  fracture  energies 
[11,12]. 

Size  effects  also  occur  at  the  macro-mechanical  level. 
For  example,  it  is  shown  in  [13]  that  the  strength  of 
notched  quasi-isotropic  composite  laminates  decreases  for 
increasing  notch  sizes  when  thin  plies  are  used.  This  effect, 
usually  known  as  the  “hole  size  effect”,  is  caused  by  the 
development  and  propagation  of  non-critical  ply-level 
damage  mechanisms  that  occur  in  the  vicinity  of  the  hole 
before  the  final  collapse  of  the  laminate.  The  exact  nature 
of  the  non-critical  damage  mechanisms  has  been  reported 
by  several  authors.  Using  Moire  interferometry  in  notched 
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[0/+45/90/— 45]s  laminates,  Mollenhauer  et  al.  [14] 
observed  a  strain  redistribution  as  a  result  of  matrix-liber 
splitting  in  the  0°  surface  ply  and  sub-surface  ply  cracking. 
Green  et  al.  [4]  reported  fiber  splitting  in  the  0°  plies,  matrix 
cracking  in  the  off-axis  plies,  and  delamination  in 
90m/— 45m/0i„]ns  carbon-epoxy  laminates  with  a  central  cir¬ 
cular  hole. 

The  observed  ply-level  damage  mechanisms  can  be 
regarded  as  a  fracture  process  zone  that  develops  before 
final  failure  of  the  laminate.  For  very  small  specimens, 
the  fracture  process  zone  affects  the  entire  width  of  the  lam¬ 
inate.  On  the  other  hand,  the  size  of  the  fracture  process 
zone  in  large  specimens  is  negligible  when  compared  with 
the  characteristic  dimensions  of  the  specimen.  The  relative 
dimension  of  the  fracture  process  zone  with  respect  to  the 
specimen  size  justifies  the  different  strengths  observed  in 
small  and  large  specimens.  Therefore,  to  predict  the  hole 
size  effect  in  quasi-brittle  materials  with  general  dimen¬ 
sions,  methods  that  account  for  the  energy  dissipated  by 
the  propagation  of  non-critical  damage  mechanisms  are 
required  [15]. 

While  the  strength  of  notched  multidirectional  laminates 
manufactured  using  thin  plies  generally  decreases  with  hole 
diameter,  Green  et  al.  [4]  reported  an  opposite  trend  for 
laminates  with  plies  with  the  same  fiber  orientation  blocked 
together  (ply-level  scaling):  for  a  4  mm  thick  [454/9O4/ 
— 454/04]s  carbon  fiber  reinforced  plastic  (CFRP)  laminate, 
increasing  the  hole  size  from  3.2  mm  to  25.4  mm  increased 
the  strength  by  51%.  This  new  finding  was  attributed  to  the 
formation  of  delaminations  at  the  edge  of  the  hole  [4].  Ply- 
blocked  specimens  exhibit  a  delamination  type  of  failure, 
and  for  small  hole  diameters  the  size  of  the  delamination 
is  relatively  large  and  grows  unstably. 

Green  et  al.  also  performed  tests  on  thickness-scaled 
CFRP  laminates  [4].  A  decrease  of  the  ultimate  strength 
with  test  specimen  thickness  was  reported  for  both  ply-level 
and  sublaminate-level  scaled  laminates,  where  the  laminate 
thickness  is  increased  by  increasing  the  number  of  sublami¬ 
nates  while  keeping  the  ply  thickness  constant.  When 
increasing  the  thickness  from  1  mm  to  8  mm,  strength 
reductions  of  16.5%  and  64.4%  were  measured  for  the 
sublaminate  level  and  ply-level  scaled  specimens,  respec¬ 
tively.  The  strength  reduction  was  attributed  to  the  higher 
energy  release  rate  at  the  interfaces  of  the  ply-level  scaled 
specimens,  which  promotes  delamination,  and  to  the  higher 
stress  concentration  relief  that  occurs  as  a  result  of  damage 
in  the  surface  plies  of  sublaminate-level  scaled  specimens. 

The  calculation  of  macro-mechanical  size  elfects  is  often 
based  on  semi-empirical  methods  that  require  calibration 
such  as  the  point  stress  and  average  stress  models  proposed 
by  Whitney  and  Nuismer  [16].  The  point  stress  model 
assumes  that  final  failure  occurs  when  the  stress  at  a  char¬ 
acteristic  distance  from  the  notch  reaches  the  unnotched 
strength  of  the  laminate.  In  the  average  stress  model,  it  is 
assumed  that  final  failure  occurs  when  the  laminate  stress 
averaged  over  a  characteristic  distance  is  equal  to  the 
unnotched  strength  of  the  laminate.  Modifications  of  the 


point  stress  and  average  stress  models  using  ply  strengths 
have  been  proposed  to  predict  the  strength  of  laminates 
with  open  and  loaded  holes  [17,18].  The  advantage  of  using 
ply  properties  rather  than  laminate  properties  is  that  the 
need  to  measure  laminate  strengths  for  every  layup  is 
avoided.  However,  the  measurement  of  the  characteristic 
distances  is  still  required  for  each  lay-up  and  geometry  [18]. 

On  the  other  hand,  models  based  on  continuum  damage 
mechanics  do  not  require  calibration,  so  they  potentially 
provide  the  means  for  a  truly  predictive  methodology  for 
the  strength  prediction  of  composite  laminates.  Continuum 
damage  models  are  defined  in  the  framework  of  the  ther¬ 
modynamics  of  irreversible  processes.  Generally  speaking, 
the  formulation  of  continuum  damage  models  starts  by 
the  definition  of  a  potential  (e.g.  the  complementary  free 
energy)  as  a  function  of  one  or  more  damage  variables  that 
is  the  basis  for  establishing  the  relation  between  the  stress 
and  the  strain  tensors.  It  is  also  required  to  define  the  dam¬ 
age  activation  functions,  i.e.  the  conditions  that  lead  to  the 
onset  of  inelastic  response,  and  the  damage  evolution  func¬ 
tions.  Some  of  the  models  proposed  in  the  literature  are 
exclusively  based  on  thermodynamic  restrictions  of  the 
constitutive  model  and  on  some  adjusting  functions  for 
damage  onset  and  evolution.  Other  models,  besides  satisfy¬ 
ing  the  thermodynamic  restrictions,  are  based  on  the  fail¬ 
ure  mechanisms  [19],  i.e.  the  damage  activation  functions 
are  related  to  the  physics  of  the  different  failure  mecha¬ 
nisms,  and  the  damage  variables  are  related  to  the  orienta¬ 
tion  of  the  ply  failure  planes  experimentally  observed. 
Mechanism-based  continuum  damage  models  can  predict 
damage  onset  and  the  extent  and  type  of  non-critical  dam¬ 
age  mechanisms.  Furthermore,  continuum  damage  models 
that  relate  the  damage  variables  to  the  normal  components 
of  the  stress  tensor  are  able  to  simulate  the  effect  of  crack 
closure  under  load  reversal  cycles.  Therefore,  such  models 
can  be  used  to  predict  the  strength  under  non-monotonic 
loading  including  load  reversals. 

The  objective  of  this  paper  is  to  investigate  the  use  of  a 
continuum  damage  model  for  the  prediction  of  size  elfects 
in  notched  carbon-epoxy  laminates  loaded  in  tension.  An 
experimental  program  is  conducted  to  measure  the  relevant 
material  properties  and  to  identify  size  effects  occurring  in 
laminates  with  different  hole  sizes.  The  recently  proposed 
continuum  damage  model  is  described  and  analysis  of  open 
hole  specimens  subjected  to  tension  loads  are  presented. 
The  analyses  results  are  compared  with  the  experimental 
data  and  with  predictions  obtained  using  a  strength  of 
materials  approach,  Linear  Elastic  Fracture  Mechanics, 
and  the  point  stress  model. 

2.  Experimental  program 

2.1.  Material  selection  and  characterization 

The  material  selected  for  the  present  study  is  Hexcel’s 
IM7-8552  carbon  epoxy  unidirectional  tape  with  a  nominal 
ply  thickness  of  0.131  mm.  The  material  was  cured  accord- 
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ing  to  the  manufacturer’s  specifications,  with  temperature 
stages  of  1 10  °C  for  1  h,  followed  by  180  °C  for  2  h.  A  pres¬ 
sure  of  7  bar  was  applied  for  the  duration  of  the  cure  cycle. 

The  fiber  volume  fraction  was  measured  using  image 
processing  techniques  resulting  in  an  average  value  of 
59.1%.  The  coefficients  of  thermal  expansion  were  mea¬ 
sured  using  a  dilatometer  and  the  resulting  values  are 
ail  =  — 5.5  X  10“V°C  for  the  longitudinal  direction,  and 
0L22  —  25.8  X  10“V°C  for  the  transverse  direction.  The  elas¬ 
tic  properties  and  strengths  were  measured  using  ASTM 
test  standards  [20-22].  Five  specimens  were  used  for  each 
test  performed. 

The  mean  measured  values  of  the  ply  elastic  properties 
are  shown  in  Table  1.  E\  and  E2  are  the  longitudinal  and 
transverse  Young’s  modulus  respectively,  G12  is  the  shear 
modulus,  and  D12  is  the  major  Poisson’s  ratio.  Table  1  also 
presents  the  standard  used  in  each  test,  the  standard  devi¬ 
ation  (STDV),  and  the  coefficient  of  variation  (CV). 

The  measured  ply  strengths  are  shown  in  Table  2.  At 
and  are  the  longitudinal  and  transverse  tensile 
strengths,  respectively.  Xc  and  Yq  are  the  longitudinal 
and  transverse  compressive  strengths,  respectively.  is 
the  in-plane  shear  strength. 

The  values  of  the  transverse  tensile  strength  ( FIJ.'*)  and  of 
the  in-plane  shear  strength  [S'^)  measured  in  the  test  spec¬ 
imens  correspond  to  the  strengths  of  unconstrained  unidi¬ 
rectional  plies.  The  transverse  tensile  and  shear  strengths 
of  constrained  plies  (in  situ  strengths)  are  higher  than  the 
ones  of  an  unidirectional  ply  [10]  and  decrease  when 
increasing  the  ply  thickness.  The  in  situ  strengths  are  calcu¬ 
lated  using  models  previously  proposed  by  the  authors, 
which  are  based  on  the  mode  I  fracture  toughness,  G2+, 
and  on  the  mode  II  fracture  toughness,  Ge  [12].  These  mod¬ 
els  use  the  simplifying  assumption  that  the  in  situ  strengths 
are  not  a  function  of  the  elastic  properties  and  geometry  of 
the  neighboring  layers. 

To  measure  the  components  of  the  fracture  toughness, 
double  cantilever  beam  (DCB)  [23]  and  four-point  bending 
end  notched  flexure  (4-ENF)  [24]  tests  were  performed.  The 
measured  components  of  the  fracture  toughness  are  shown 
in  Table  3. 

The  in  situ  strengths  are  calculated  as  functions  of  the 
fracture  toughness  and  ply  elastic  properties  using  the 
models  described  in  [12]  with  a  shear  response  factor 
fi  =  2.98  X  10“*  MPa“^.  The  calculated  in  situ  strengths 
are  shown  in  Table  4. 

The  shear  strength  in  the  transverse  direction  is  calcu¬ 
lated  as  [25,26] 


Table  1 

Measured  ply  elastic  properties  for  IM7-8552 


Property 

Standard 

Mean  value 

STDV 

CV  (%) 

(GPa) 

Ref.  [20] 

171.42 

2.38 

1.39 

E2  (GPa) 

Ref.  [20] 

9.08 

0.09 

1.03 

G,2  (GPa) 

Ref.  [22] 

5.29 

0.13 

2.53 

«|2 

Ref.  [20] 

0.32 

0.02 

6.18 

Table  2 


Measured  ply  strengths  for  IM7-8552 

Property 

Standard 

Mean  value  (MPa) 

STDV  (MPa) 

CV  (%) 

TV 

Ref  [20] 

2326.2 

134.1 

5.8 

Ref  [21] 

1200.1 

145.7 

12.1 

vud 

Ij 

Ref  [20] 

62.3 

5.3 

8.5 

Yc 

Ref  [21] 

199.8 

20.5 

10.2 

Sl^ 

Ref  [22] 

92.3 

0.6 

0.7 

Table  3 

Measured  fracture  energies  for  transverse  fracture  for  IM7-8552  (kj/m^) 

Property 

Mean  value 

STDV 

CV  (%) 

G2+ 

0.2774 

0.0246 

0.88 

Ge 

0.7879 

0.0803 

10.19 

Table  4 

Calculated  in  situ  strengths  for  IM7-8552  (MPa) 

Ply  configuration 

5V 

Thin  embedded  ply 

160.2 

130.2 

Thin  outer  ply 

101.4 

107.0 

„  T,  7  •  COSao  \  ,,, 

=  Fccosao  smao +;^ — ^  (1) 

where  ao  is  the  fracture  angle  of  a  ply  under  pure  transverse 
compression  [27].  For  a  fracture  angle  ag  =  53°,  the  shear 
strength  in  the  transverse  direction  is  calculated  as 
5t  =  75.3  MPa. 

The  continuum  damage  model  also  requires  the  fracture 
energies  per  unit  surface  for  longitudinal  failure,  Gi+  (ten¬ 
sion)  and  Gi_  (compression).  These  energies  were  mea¬ 
sured  using  the  Compact  Tension  (CT)  and  Compact 
Compression  (CC)  tests  in  cross-ply  laminates  proposed 
by  Pinho  et  al.  [28,29].  The  measured  fracture  energies 
per  unit  surface  are  shown  in  Table  5. 

2.2.  Notched  laminates 

Tests  of  notched  composite  laminates  were  performed  to 
quantify  the  size  effect  and  to  obtain  empirical  data  to  val¬ 
idate  the  numerical  model.  Quasi-isotropic  laminates  were 
manufactured  in  Hexcel  IM7-8552  CFRP  with  a  stacking 
sequence  of  [90/0/±45]3s. 

The  unnotched  tensile  strength  of  the  laminate,  X\,  was 
measured  using  five  test  specimens  and  the  average  value 
obtained  was  845.1  MPa.  The  average  value  of  the  failure 
strain,  12,900;U8,  was  measured  in  the  five  test  specimens 
using  strain  gages. 


Table  5 

Measured  fracture  energies  for  longitudinal  fracture  for  IM7-8552  (kJ/m^) 
Property  Mean  value  STDV  CV  (%) 

Gi+  81.5  6.1  7.6 

Gi_  106.3  2.2  2.1 
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The  notched  test  specimens  were  machined  using  a  pro¬ 
cedure  that  prevents  delaminations  in  the  regions  close  to 
the  insertion  point  and  the  exit  of  the  drill  bit.  Sacrificial 
frontal  and  backing  plates  were  used  to  clamp  the  speci¬ 
mens  during  the  drilling  process.  All  test  specimens  were 
machined  to  class  1  hole  quality  used  in  aerospace  [30]. 
No  damage  was  observed  in  a  sample  of  test  specimens 
inspected  using  X-rays. 

Specimens  with  five  different  hole  diameters,  d=2  mm, 
4  mm,  6  mm,  8  mm,  10  mm  and  with  a  width-to-diameter 
ratio  wid  equal  to  6  were  tested  in  a  MTS  servo-hydraulic 
machine  following  the  ASTM  D-5766  standard  [31] 
according  to  the  test  matrix  shown  in  Table  6.  Five  speci¬ 
mens  were  tested  for  each  geometry. 

The  specimens  labeled  OHT3,  OFIT6  and  OHT9  were 
instrumented  with  two  strain  gages  in  the  positions  sche¬ 
matically  shown  in  Fig.  1. 

The  distances  d^  shown  in  Fig.  1  are  respectively 
10.5  mm,  13.5  mm,  and  12.5  mm  for  the  test  specimens 
OHT3,  OHT6  and  OFIT9.  The  specimens  OHTIO  and 
OHTll  were  not  instrumented.  Acoustic  emission  (AE) 
sensors  were  used  in  one  test  specimen  for  each  size. 


Table  6 

Open  hole  tension  test  matrix 


Specimen  ref. 

d  (mm) 

w  (mm) 

w/d 

OHTll 

2 

12 

6 

OHTIO 

4 

24 

6 

OHT3 

6 

36 

6 

OHT6 

8 

48 

6 

OHT9 

10 

60 

6 

Fig.  2  shows  the  applied  load  and  the  cumulative  num¬ 
ber  of  AE  signals  as  a  function  of  time  for  one  OHT3  test 
specimen. 

From  the  AE  signals  shown  in  Fig.  2,  it  can  be  con¬ 
cluded  that  non-critical  damage  mechanisms  accumulate 
well  before  final  failure  of  the  specimen,  creating  a  fracture 
process  zone  (FPZ).  Similar  results  are  observed  in  the 
OHT6  and  OFIT9  specimens,  as  well  as  in  other  experimen¬ 
tal  investigations  [4,14]. 

The  remote  failure  stress  is  defined  using  the  failure  load 
measured  in  the  tests  {P)  and  the  measured  values  of  the 
specimen  thickness  (A)  and  width  (w)  as:  The 

remote  failure  stresses  obtained  for  the  different  geometries 
are  summarized  in  Table  7. 

The  failure  mode  observed  in  all  specimens  is  net-section 
tension,  as  shown  in  Fig.  3.  Fig.  4  shows  the  relation 
between  the  remote  stress  and  the  strain  measured  by  strain 
gages  SG3  for  one  test  specimen  of  each  of  the  three  differ¬ 
ent  geometries. 

The  experimental  results  presented  in  Table  7  clearly 
identify  a  size  effect;  an  increase  in  the  hole  diameter  from 
2  mm  to  10  mm  results  in  a  32.8%  reduction  in  the  strength. 
The  observed  size  effect  is  caused  by  the  development  of  the 
fracture  process  zone  identified  in  the  AE  results,  which  re¬ 
distributes  the  stresses  and  dissipates  energy.  In  small  spec¬ 
imens,  the  fracture  process  zone  extends  towards  the  edges 
of  the  specimen  and  the  average  stress  at  the  fracture  plane 
tends  to  the  unnotched  strength  of  the  laminate. 


Fig.  2.  Applied  load  and  AE  signals  as  a  function  of  time  for  the  specimen 
with  a  6  mm  diameter  hole. 


Table  7 


Results  of  open-hole  tensile  tests 

Hole  diameter  (mm) 

(MPa) 

STDV  (MPa) 

CV  (%) 

2 

555.7 

15.3 

2.8 

4 

480.6 

21.4 

4.5 

6 

438.7 

25.3 

5.8 

8 

375.7 

15.1 

4.0 

10 

373.7 

14.1 

3.8 

Fig.  1 .  Position  of  strain  gages. 
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Fig.  3.  Net-section  tension  failures  in  specimens  with  wld=  6. 


o“,MPa 


Fig.  4.  Relation  between  remote  stress  and  longitudinal  strain  in  SG3. 


The  effect  of  size  on  the  strength  can  be  explained  using 
a  simple  example  based  on  the  cohesive  crack  model,  which 
is  well-suited  to  simulate  fracture  of  quasi-brittle  materials 
[32].  Consider  that  the  fracture  process  zone  is  represented 
by  a  cohesive  crack  with  the  simple  constitutive  relation 
shown  in  Fig.  5a. 

The  cohesive  constitutive  model  relates  the  laminate 
cohesive  stress,  a,  to  the  crack  opening,  w,  and  must  satisfy 
the  following  condition:  f^(T(w)dw=  Gc-  Structural  col¬ 
lapse  occurs  when  a  point  along  the  fracture  plane  reaches 
the  critical  opening,  w^,  and  the  corresponding  length  of 


a 

a 


b 


tnttttttttttmt 


w 


— > 

c  W 
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mmi 


mim 


Fig.  5.  Cohesive  crack  constitutive  law  and  fracture  process  zone,  (a) 
Constitutive  model  and  (b)  stress  distributions  along  the  fracture  plane. 


the  fracture  process  zone  can  be  estimated  using  the  Irwin 
model  as  /fpz  «  ^  [33,34]. 

Based  on  the  constitutive  law  shown  in  Fig.  5a,  it  is  pos¬ 
sible  to  schematically  represent  the  stress  distribution  at 
failure  along  the  fracture  planes  of  specimens  with  different 
sizes,  as  shown  in  Fig.  5b.  It  is  observed  that  in  small  spec¬ 
imens  the  fracture  process  zone  extends  towards  the  edges, 
whereas  in  large  specimens  the  fracture  process  zone  is  con¬ 
fined  to  the  vicinity  of  the  hole.  As  a  consequence,  the  aver¬ 
age  stress  acting  on  the  fracture  plane,  and  hence  the 
strength,  are  larger  for  small  specimens. 


3.  Simulation  of  the  effect  of  size  on  strength 


Strength  prediction  methods  uniquely  based  on  stress  or 
strain  failure  criteria  are  unable  to  predict  the  size  effects 
observed  in  notched  specimens.  Consider  for  example  a  cal¬ 
culation  of  the  final  failure  of  a  specimen  with  a  central 
hole  using  the  value  of  the  longitudinal  stress  in  the  fiber 
direction  (maximum  stress  criterion).  The  distribution  of 
the  longitudinal  stress  in  the  critical  plies,  the  0°  plies  along 
the  fracture  plane,  defined  by  0  =  90°  in  Fig.  1,  can  be  cal¬ 
culated  using  an  approximate  closed-form  solution  as  [35] 

(Til  =  (7xx{d,y){Qna\^  +  Q\ia\2)  (2) 


where  a*-  are  the  components  of  the  laminate  compliance 
matrix  defined  as  [36] 


[a*]  =  ti\A]  ‘ 


(3) 


where  the  matrix  [A]  relates  the  in-plane  forces  per  unit 
length  to  the  mid-plane  strains.  are  the  components  of 
the  plane  stress  transformed  reduced  stiffness  matrix  of 
the  0°  plies  [18],  and  A  is  the  thickness  of  the  laminate. 

The  through-the-thickness  averaged  normal  stress  in  the 
fracture  plane  for  a  quasi-isotropic  laminate  is  calculated 
by  Tan  [35]  as 


y  ^  d/2 


2  +  (1  —  djw) 
6(1  —  d/w) 


(4) 


where  is  the  remote  tensile  stress. 

From  Eqs.  (2)  and  (4)  it  is  clear  that  for  the  same  mate¬ 
rial  and  stacking  sequence  the  stress  concentration  factor. 
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and  hence  the  maximum  longitudinal  stress  in  the  0°  ply, 
depends  on  the  ratio  between  the  specimen  hole  diameter 
and  width.  Applying  the  maximum  stress  criterion  and 
using  Eqs.  (2)  and  (4) 

-  ,5) 

At  [2  + (\  -  d/w)  ]{Qna\^+ Qna\2) 

Eq.  (5)  demonstrates  that  the  application  of  the  maximum 
stress  criterion  results  in  the  same  strength  prediction  for 
different  hole  diameters  when  the  d/w  ratio  is  held  constant. 
The  lack  of  size  effect  on  the  predicted  strength  clearly  con¬ 
tradicts  the  experimental  observations. 

3.1.  Linear  Elastic  Fracture  Mechanics 


account  that  the  hole  radius  tends  to  zero,  in  which  case 
the  function  f(a,  R)  tends  to  one,  leaving 

(10) 

where  Xj  is  the  tensile  strength  of  the  unnotched  laminate. 

From  (8)  and  (10),  the  equation  proposed  by  Waddoups 
et  al.  [13]  is  obtained 

a^=X\/f{a,R)  (11) 

The  strength  of  the  laminate  containing  an  open-hole  is  pre¬ 
dicted  using  two  parameters:  the  length  of  the  inherent  flaw, 
a,  that  needs  to  be  calculated  from  a  baseline  specimen,  and 
the  unnotched  tensile  strength  of  the  laminate,  Xj. 

3.2.  Point-stress  model 


There  are  two  approaches  that  can  be  used  with  Linear 
Elastic  Fracture  Mechanics  (LEFM)  to  calculate  the  effect 
of  size  on  the  strength  of  notched  composite  laminates.  In 
the  first  approach,  it  is  assumed  that  the  length  u  of  a  pre¬ 
existing  crack  in  the  laminate  is  scaled  in  the  same  propor¬ 
tion  of  the  hole  diameter  and  specimen  width  and  that  the 
critical  value  of  the  laminate’s  stress  intensity  factor,  Ki^  is 
independent  of  the  crack  length.  Consider  two  specimens 
with  hole  diameters  d^  and  ^2-  The  stress  intensity  factor 
at  failure  is 


Ai,  =  afF 


Taking  into  account  the  fact  that  the  crack  length  is  pro¬ 
portional  to  the  hole  diameter  and  that  the  finite  width  cor¬ 
rection  factors,  F{wld,ald),  are  equal  for  scaled  geometries, 
the  failure  stress  of  a  specimen  with  a  hole  diameter  ^2  can 
be  calculated  from  the  failure  stress  of  the  specimen  with  a 
hole  diameter  di 


The  point-stress  model  (PSM)  proposed  by  Whitney  and 
Nuismer  [16],  considers  that  ultimate  failure  occurs  when 
the  stress  at  a  given  distance  from  the  hole  boundary,  t'o,, 
reaches  the  unnotched  strength  of  the  laminate,  Xj.  An 
alternative  version  of  the  point  stress  model  uses  the  ply 
stresses  and  strengths,  so  that  it  is  not  necessary  to  measure 
the  strength  for  every  different  laminate. 

Using  Eqs.  (2)  and  (4),  the  strength  predicted  using  the 
PSM  is 


=At 


^2+(i-^r 

,  6(1 


d  +  Ivo 


+  3 


d  +  Ivo 


X  (fill's'll  +  612012) 


(12) 


Failure  is  predicted  using  two  parameters:  the  characteris¬ 
tic  distance  in  tension  and  the  longitudinal  tensile 
strength  of  the  ply,  X^. 


(V) 

The  second  approach  to  predict  size  effects  using  LEFM  is 
the  inherent  ffaw  model  (IFM)  proposed  by  Waddoups 
et  al.  [13].  It  is  considered  that  the  non-critical  damage 
mechanisms  occurring  before  ultimate  failure  of  a  compos¬ 
ite  laminate  can  be  lumped  into  a  constant  “region  of  in¬ 
tense  energy”,  or  “inherent  ffaw”,  of  length  a.  The 
critical  value  of  the  stress  intensity  factor  of  a  plate  with 
a  hole  of  radius  R  is  given  by 

K\^,=  f{a,R)d°°y/na  (8) 

where  f(a,R)  is  Bowie’s  solution  for  the  calculation  of  the 
stress  intensity  factor  of  two  cracks  emanating  from  a  cir¬ 
cular  hole,  given  as  [37,38] 


f{a,R)  =  0.5 


a  \ 
d/2-\-  a) 


1  1.243 


a  \ 
d/2  +  a) 


(9) 


Waddoups  et  al.  [13]  considered  that  the  strength  of  an 
unnotched  specimen  can  be  predicted  by  taking  into 


3.3.  Continuum  damage  model 

Continuum  damage  mechanics  is  a  methodology  well 
suited  for  the  simulation  of  damage  evolution  and  ultimate 
failure  of  composites  under  general  loads  and  boundary 
conditions  for  which  no  analytical  solution  is  available. 
The  continuum  damage  model  used  here  is  based  on  previ¬ 
ous  work  by  the  authors  [19,39,40].  The  main  aspects  of  the 
continuum  damage  model  are  presented  in  the  following 
sections.  The  full  details  of  the  model  can  be  found  in  Refs. 
[19,39,40]. 

3.3.1.  Constitutive  model 

The  proposed  definition  for  the  complementary  free 
energy  density  of  a  ply  is 

^  <^11  I  a2  <^12 

2(l-di)Ci  2(l-d2)C2  Cl  "  ''■^2(l-d6)Gi2 
+  (ofiiO'ii  +  ina22)S.T  H-  (^nffn  -t-  (13) 

where  the  damage  variable  d  1  is  associated  with  longitudi¬ 
nal  (fiber)  failure,  d2  is  the  damage  variable  associated  with 
transverse  matrix  cracking,  and  ds  is  the  damage  variable 
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associated  with  longitudinal  and  transverse  cracks.  /In  and 
^22  are  the  coefficients  of  hygroscopic  expansion  in  the  lon¬ 
gitudinal  and  transverse  directions,  respectively.  ST  and 
AM  are  the  differences  of  temperature  and  moisture  con¬ 
tent  with  respect  to  the  corresponding  reference  values. 
The  coefficients  of  thermal  expansion  of  a  ply  are  also  af¬ 
fected  by  the  failure  mechanisms.  The  exact  dependence 
of  the  coefficients  of  thermal  expansion  with  damage  can 
be  obtained  for  simple  laminates  in  the  absence  of  stress 
gradients  [41].  These  conditions  are  not  met  by  the  lami¬ 
nate  under  investigation  here  and  the  effects  of  damage 
on  the  coefficients  of  thermal  expansion  are  neglected. 

The  strain  tensor  is  equal  to  the  derivative  of  the  com¬ 
plementary  free  energy  density  with  respect  to  the  stress 
tensor 


£ 


W 

0(7 


H  ;  (7  -1-  aST  +  pAM 


(14) 


The  lamina  compliance  tensor  can  be  represented  as 


1  _  in2  n 

(l-di)£i  £,  " 

_  H12.  1  n 

El  (l-d2)E2  " 
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0 


1 

(l-d6)Gi2  _ 


(15) 


The  closure  of  transverse  cracks  under  load  reversal  is  ta¬ 
ken  into  account  by  defining  four  damage  variables  associ¬ 
ated  with  longitudinal  and  transverse  damage.  To 
distinguish  between  the  active  and  the  passive  damage  vari¬ 
ables,  it  is  necessary  to  define  the  longitudinal  and  trans¬ 
verse  damage  modes  as  follows: 


di  =  di+ 

d2  =  d2+ 


(gii) 

kill 

(0'22) 

k22| 


+  di_ 

+  di2- 


(-gii) 

kill 

{—(J22) 

k22| 


(16) 


where  (x)  is  the  McCauley  operator  defined  as 
(x)  :=  (x  +  |x|)/2. 


3.3.2.  Damage  activation  functions 

The  determination  of  the  domain  of  elastic  response 
under  complex  stress  states  is  an  essential  component  of 
an  accurate  damage  model.  It  is  assumed  that  the  elastic 
domain  is  enclosed  by  four  surfaces,  each  of  them  account¬ 
ing  for  one  damage  mechanism:  longitudinal  and  trans¬ 
verse  fracture  under  tension  and  compression.  Those 
surfaces  are  formulated  by  the  damage  activation  functions 
based  on  the  LaRC04  failure  criteria  [26]. 

The  four  damage  activation  functions,  T’at,  associated 
with  damage  in  the  longitudinal  {N  =  1+,  1—)  and  trans¬ 
verse  {N  =  2+,  2—)  directions  represented  in  Fig.  6,  are 
defined  as 


Fi+ =  -  ri+ ^  0;  -ri_  <  0 

^2+  =  (/'2+  -  ^2+  <  0;  F2-  =  (/)2-  -  7-2-  <  0 


(17) 


where  the  loading  functions  (N  —  1+,  1—,  2+,  2—)  de¬ 
pend  on  the  strain  tensor  and  material  constants  (elastic 
and  strength  properties).  The  elastic  domain  thresholds 

{N  =  1+,  1—,  2+,  2—)  take  an  initial  value  of  1  when 
the  material  is  undamaged,  and  they  increase  with  damage. 
The  elastic  domain  thresholds  are  related  to  the  damage 
variables  (M=  1+,  1—,  2+,  2—,  6)  by  the  damage  evo¬ 
lution  laws. 

The  current  values  of  the  elastic  domain  thresholds 
are  obtained  using  the  loading  functions  fn/  according  to 
the  following  equations  [19,39,40]: 


ri+  =  max 1,  max  {(/)]^},  max  {(/)']_} 
ri_  =  max  1 1,  max  {(/)]  }  I 
r2+  =  max  1 1 ,  max  { (/)2_  } ,  max{  (f)2+  }  | 
r2_  =  max  1 1,  max  {k2-}| 


(18) 


Fig.  6.  Fracture  surfaces  and  corresponding  internal  variables,  (a)  Longitudinal  tensile  fracture,  (b)  longitudinal  compressive  fracture,  (c)  transverse 
fracture  with  oc  =  0°  and  (d)  transverse  fracture  with  a.  =  53°. 
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3.3.2. 1.  Longitudinal  tensile  fracture.  The  LaRC04  crite¬ 
rion  for  fiber  tension  is  defined  as 


01+  = 


(Til  —  VnOii 
Xj 


(19) 


where  the  effective  stress  tensor  d  is  computed  as 
0-  =  Ho  ‘  :  £.  Ho  is  the  undamaged  compliance  tensor. 


3. 3. 2. 2.  Longitudinal  compressive  fracture.  The  damage 
activation  function  used  to  predict  damage  under  longitu¬ 
dinal  compression  (un  <  0)  and  in-plane  shear  (fiber  kink¬ 
ing)  is  established  as  a  function  of  the  components  of  the 
stress  tensor  in  a  coordinate  system  [m)  representing 
the  fiber  misalignment 


01- 


(20) 


where  the  coefficient  of  longitudinal  influence  can  be 
approximated  as  [26] 


5l  cos(2(3Co) 
Yc  COS^  OLq 


(21) 


with  Ko  =  53°  [27],  The  components  of  the  effective  stress 
tensor  in  the  coordinate  system  associated  with  the  rotation 
of  the  fibers  are  calculated  as 


^22  =  ^11  +  ^22  cos^  (p^  —  2\d\2  \  sin  cos  qf 

^12  =  (^22  —  ffii)  sin  qf  cos  qf  +  lu^Kcos^  qf  —  sin^  q)'^) 

(22) 


The  misalignment  angle  [qf)  is  determined  using  standard 
shear  and  longitudinal  compression  strengths,  Sl  and  Xc, 
respectively  [26] 


3. 3. 2. 3.  Transverse  fracture  perpendicular  to  the  mid-plane 
of  the  ply.  Transverse  matrix  cracks  perpendicular  to 
the  mid-plane  of  the  ply,  i.e.  with  ao  =  0°,  are  created  by 
a  combination  of  in-plane  shear  stresses  and  transverse  ten¬ 
sile  stresses,  or  in-plane  shear  stresses  and  small  transverse 
compressive  stresses.  These  conditions  are  represented  by 
the  following  failure  criteria: 


where  g  is  the  fracture  toughness  ratio  defined  as  g  = 


3. 3. 2. 4.  Transverse  compressive  fracture.  The  matrix  fail¬ 
ure  criterion  for  transverse  compressive  stresses  consists 
of  a  quadratic  interaction  between  the  effective  shear  stres¬ 
ses  acting  on  the  fracture  plane 


02- 


if  $22  <  0 


(25) 


where  the  effective  stresses  f  and  are  computed 
^ff  =  (-ff22cos(ao)(sin(ao)  -  J7’^cos(ao)  cos(0s))) 
T^ff  =  (cos(ao)(|5i2|  +  f;^ff22  cos((Xo)  sin(0J)) 


with  if  = 


-1 


tan(2ao) 


and  0,  =  arctan 


(  -bi2l  ^ 
\^S22  sin{txi,)  J  ■ 


as  [26] 
(26) 


3.3.3.  Damage  evolution  laws  and  numerical  implementation 

Strain-softening  constitutive  models  that  do  not  take 
into  account  the  finite  element  discretization  produce 
results  that  are  mesh-dependent,  i.e.  the  solution  is  non¬ 
objective  with  respect  to  the  mesh  refinement  and  the  com¬ 
puted  energy  dissipated  decreases  with  a  reduction  of  the 
element  size  [42,43].  An  effective  solution  to  assure  objec¬ 
tive  solutions  consists  of  using  a  characteristic  length  of 
the  finite  elements  (f)  in  the  definition  of  the  constitutive 
model  [42].  As  schematically  shown  in  Fig.  7,  the  post-peak 
response  of  the  material  is  scaled  as  a  function  of  the  ele¬ 
ment  size  to  keep  the  computed  energy  dissipation  indepen¬ 
dent  of  the  size  of  the  element,  and  equal  to  the  material 
fracture  energy. 

The  energetic  regularization  of  the  model  proposed 
requires  the  fracture  energies  associated  with  the  four  frac¬ 
ture  planes  shown  in  Fig.  6.  These  fracture  energies  were 
measured  in  the  experimental  program  and  are  used  in 
the  damage  evolution  laws. 

The  exponential  damage  evolution  laws  proposed  by  the 
authors  [19,39,44]  are  expressed  in  the  following  general 
form: 

dw  =  1  -  .  /  ,  exp{^;,^[l  -fN{rN)\}f{rK)  (27) 

Jnvn) 

where  the  function  fNfN)  is  selected  to  force  the  softening 
of  the  constitutive  relation  and  it  is  taken  as  being  indepen¬ 
dent  of  the  material.  The  term/jr^:)  represents  the  coupling 
factor  between  damage  laws  and  elastic  threshold  domains. 
The  specific  damage  evolution  laws  for  each  damage  vari¬ 
able  are  presented  in  [19,39^4]. 

The  regularization  of  the  energy  dissipated  is  performed 
by  integrating  the  rate  of  energy  dissipation  for  each  failure 


Fig.  7.  Scaling  of  constitutive  model  for  different  element  sizes. 
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mode.  The  energy  dissipated  in  each  failure  mode  must  be 
independent  of  the  element  size,  and  must  be  equal  to  the 
fracture  energy  measured  in  the  experiments 


r 

0dM  dru  r  ’ 


M  —  1+,  1  — ,  2+,  2 — ,  6 


(28) 


Using  (27)  in  (28),  it  is  possible  to  numerically  integrate  the 
resulting  equation  and  calculate  the  parameters  A  m  that  as¬ 
sure  a  mesh-independent  solution  [19], 

The  constitutive  model  was  implemented  in  the  ABA- 
QUS  Finite  Element  (FE)  code  [45]  as  a  user-written 
UMAT  subroutine. 


33.4.  Mesh  objectivity  and  unidirectional  notched  specimen 

The  mesh  objectivity  of  the  model  proposed  is  illustrated 
by  simulating  the  response  of  a  notched  [90]24  CFRP  lam¬ 
inate  loaded  in  tension.  The  specimen  simulated  is  150  mm 
long,  12  mm  wide,  3  mm  thick,  and  contains  a  central  cir¬ 
cular  notch  with  a  diameter  of  6  mm.  The  properties  used 
are  reported  in  Tables  1-3. 

Two  FE  models  with  different  mesh  refinements  and 
using  the  damage  model  outlined  in  the  previous  sections 
were  created.  Models  1  and  2  use,  respectively,  6  and  20  ele¬ 
ments  along  the  fracture  plane.  Only  one-half  of  the  spec¬ 
imen  width  is  modeled.  The  details  of  the  two  meshes  are 
shown  in  Fig.  8. 

Fig.  9  shows  the  load-displacement  relation  predicted 
using  the  constitutive  model  proposed.  It  is  observed  that 
the  solution  is  independent  of  the  mesh  refinement. 

In  order  to  demonstrate  the  error  introduced  by  not 
accounting  for  element  size,  two  analyses  with  different  lev¬ 
els  of  mesh  refinement  were  also  conducted  with  a  constitu¬ 
tive  model  that  is  not  adjusted  using  Eq.  (28).  Instead,  a 
constant  softening  parameter  A2+  =  1.5  is  used,  indepen¬ 
dently  of  the  mesh  refinement.  The  load-displacement  rela¬ 
tion  predicted  by  this  model  is  shown  in  Fig.  10.  It  is  clear 
from  this  figure  that  the  maximum  load  and  energy  dissipa¬ 
tion  predicted  are  a  function  of  the  refinement  of  the  mesh. 

3.3.5.  Quasi-isotropic  open  hole  tension  specimens 

Finite  element  models  of  all  OHT  specimens  shown  in 
the  test  matrix  presented  in  Section  2  were  created  using 


a 


b 


Fig.  8.  Different  mesh  refinements:  (a)  mesh  1;  (b)  mesh  2. 


Fig.  9.  Load-displacement  relation  predicted  using  the  model  proposed. 


Fig.  10.  Load-displacement  relation  predicted  using  the  non-adjusted 
model. 


ABAQUS  [45]  four-node  S4  shell  elements.  The  difference 
between  the  working  and  reference  temperatures  used  to 
calculate  the  residual  thermal  stresses  was  —  155°C.  An 
implicit  dynamic  analysis  was  subsequently  performed, 
and  the  loading  rate  used  in  the  experiments,  2  mm/min, 
was  also  applied  to  the  numerical  models.  The  use  of  an 
implicit  dynamic  finite  element  model  enables  the  predic¬ 
tion  of  the  load  drop  that  occurs  when  the  specimens  fail 
catastrophically.  The  material  properties  used  are  pre¬ 
sented  in  Tables  1-5. 

Delamination  is  not  simulated  by  the  model.  As 
explained  by  Green  at  al.  [4],  delamination  is  the  driving 
failure  mechanism  for  ply-blocked  laminates,  but  not  for 
sublaminate-level  scaled  laminates,  such  as  those  used  in 
this  work.  The  simulation  of  delamination  is  required  for 
ply-blocked  laminates,  and  can  be  performed  using  cohe¬ 
sive  elements  connecting  several  shell  elements  that  repre¬ 
sent  the  layers  [46]. 

The  models  simulate  the  fracture  process  from  the  onset 
of  damage  up  to  structural  collapse.  Fig.  1 1  shows  the  evo¬ 
lution  of  fiber  fracture  predicted  in  the  top  0°  ply,  as  well 
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Fig.  11.  Evolution  of  fiber  fracture  in  the  top  0°  ply  for  the  specimen 
OHT9. 

the  plane  of  localized  deformation  (fracture  plane)  for 
specimen  configuration  OHT9. 

Fig.  12  shows  the  relation  between  the  applied  remote 
stress  and  the  longitudinal  deformation  measured  using 
the  strain  gages  and  the  corresponding  numerical  predic¬ 
tions  in  the  specimen  OHT6.  The  numerical  results  corre- 


Fig.  12.  Experimental  and  numerical  results-  specimen  OHT6,  strain 
gages  SG2  and  SG3. 


Table  8 


Comparison  between  experimental  and  numerical  failure  stresses  (MPa) 


Flole  diameter  (mm) 

(7“,  experimental 

(7“,  numerical 

Error  (%) 

2 

555.7 

553.6 

-0.4 

4 

480.6 

463.0 

-3.7 

6 

438.7 

430.0 

-2.0 

8 

375.7 

415.0 

-t-10.5 

10 

373.7 

405.6 

+8.5 

spond  to  the  averaged  strain  calculated  in  the  group  of 
elements  whose  position  and  total  area  correspond  approx¬ 
imately  to  the  area  where  the  strain  gages  were  bonded  to 
the  specimen.  The  location  of  the  different  strain  gages  is 
shown  in  Fig.  1. 

The  remote  failure  stresses  measured  in  the  experimental 
program  and  predicted  by  the  numerical  model  are  shown 
in  Table  8. 

From  the  comparison  between  the  experimental  and 
numerical  results,  both  in  terms  of  stress-strain  relations 
and  failure  stresses,  it  can  be  concluded  that  the  model  is 
capable  of  predicting  with  good  accuracy  the  response  of 
all  OHT  specimens  that  were  tested. 

3.4.  Comparison  of  approaches 

The  four  methods  previously  described,  i.e.  strength  of 
materials,  LEFM-scaled,  LEFM-inherent  flaw  model, 
point  stress  model,  and  continuum  damage  model  were 
applied  to  predict  the  size  effect  for  the  specimens  described 
in  Section  2.2. 

Eq.  (7)  provides  the  LEFM-scaled  prediction  for  the 
notched  strength  of  the  laminate  when  all  the  in-plane 
dimensions  are  scaled.  The  average  failure  stress  measured 
in  the  specimens  with  a  hole  diameter  of  6  mm  was  used  in 
the  LEFM  model  to  predict  the  strength  of  the  specimens 
with  different  geometries. 

Eq.  (11)  provides  the  LEEM-inherent  flaw  model  pre¬ 
diction  of  the  notched  strength.  The  specimen  with  a 
6  mm  hole  diameter  is  used  to  calculate  the  length  of  the 
inherent  flaw.  Using  the  measured  mean  failure  stress  in 
Eq.  (11),  the  length  of  the  inherent  flaw  is  calculated  as 
a  —  1.28  mm. 

The  point-stress  prediction  of  the  size  effect  is  performed 
using  Eq.  (12).  The  characteristic  distance  of  0.75  mm  was 
obtained  by  using  the  measured  mean  failure  stress  in  the 
specimen  with  a  6  mm  hole  diameter.  This  value  of  the 
characteristic  distance  is  used  to  predict  the  strength  of 
the  other  specimens. 

The  predictions  of  the  normalized  strength  as  a  function 
of  the  hole  diameter  obtained  using  the  different  models  are 
shown  in  Fig.  13. 

It  can  be  observed  that  both  the  point  stress  and  LEFM- 
IFM  models  can  predict  with  reasonable  accuracy  the  size 
effect  law  of  notched  composite  laminates.  The  point  stress 
and  inherent  flaw  models  are  particularly  accurate  for  spec¬ 
imens  with  hole  diameters  close  to  the  diameter  used  to  cal¬ 
culate  the  characteristic  distance  (PSM)  and  the  length  of 
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Laminate  unnotchcd  strength, 


ZA  - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 

0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1.0  1.1 

log{d) 

Fig.  13.  Predictions  of  size  effects  in  CFRP  plates  with  wid  =  6. 

the  inherent  flaw.  For  specimens  with  small  hole  diameters, 
the  predictions  lose  accuracy.  Therefore,  to  accurately  pre¬ 
dict  the  notched  strength  of  laminates  these  models  require 
the  calculation  of  the  characteristic  distance  and  length  of 
inherent  flaw  for  different  geometries,  and  the  deflnition 
of  an  extrapolation  procedure  to  deflne  the  values  of  these 
parameters  for  other  geometries  [18].  It  should  also  be 
noted  that  the  basic  equation  used  in  the  inherent  flaw 
model,  Eq.  (11),  is  only  valid  when  finite  width  effects  are 
negligible,  which  is  the  case  of  the  specimens  tested.  For 
smaller  ratios  between  the  specimen  width  and  hole  diam¬ 
eter,  the  inherent  flaw  model  should  be  modified. 

The  continuum  damage  model  can  predict  the  size  effect 
law  observed  in  the  experiments,  especially  for  specimens 
with  hole  diameters  smaller  than  6  mm.  Unlike  the  point 
stress  and  inherent  flaw  models,  the  continuum  damage 
model  does  not  require  any  adjustment  parameter  and  only 
uses  material  properties  that  are  measured  at  the  ply  level 
as  well  as  the  fracture  energies. 

Fig.  1 3  indicates  that  the  use  of  the  LEFM-scaled  model 
results  in  accurate  predictions  for  hole  sizes  between  6  mm 
and  10  mm.  However,  the  strength  is  overpredicted  for 
small  hole  diameters.  Eor  small  specimens,  the  damaged 
region  in  the  vicinity  of  the  hole  cannot  be  considered  to 
be  negligible  when  compared  with  the  characteristic  dimen¬ 
sions  of  the  specimen,  and  LEFM  is  not  applicable. 

LEFM-scaled  predictions  are  also  inaccurate  for  large 
specimens  because  the  notched  strengths  of  those  speci¬ 
mens  tend  to  a  constant  value  [4].  Bazant  [15]  relates  this 
asymptotic  structural  response  to  the  invariance  of  the  size 
of  the  fracture  process  zone  when  the  characteristic  dimen¬ 
sions  of  large  specimens  are  increased.  It  should  also  be 
noted  that  the  LEFM  predictions  based  on  scaled  speci¬ 
mens  always  result  in  a  line  with  a  —1/2  slope  that  passes 
through  the  baseline  point  (Fig.  13).  This  means  that  the 
use  of  a  small  hole  diameter  as  the  baseline  point  would 
result  in  severe  underpredictions  of  the  notched  strength 
of  larger  specimens. 


The  maximum  stress  criterion  for  longitudinal  failure  is 
unable  to  predict  size  effects  and  always  underpredicts  the 
strength  of  notched  laminates.  For  a  hole  diameter  of 
2  mm,  the  application  of  the  maximum  stress  criterion 
results  in  an  error  of  —49.1%.  The  error  associated  with 
the  strength  of  materials  approach  is  even  larger  when 
using  a  failure  criterion  for  transverse  (matrix)  cracking, 
which  occurs  before  fiber  fracture,  or  failure  criteria  that 
are  unable  to  distinguish  fiber  and  matrix  failure  modes. 

4.  Conclusions 

The  size  effect  in  notched  IM7-8552  CFRP  was  identified 
and  quantified  in  an  experimental  program.  The  acoustic 
emission  results  show  that  final  fracture  is  preceded  by  a  pro¬ 
cess  of  accumulation  of  non-critical  damage  mechanisms. 

By  comparing  the  experimental  data  with  the  different 
models  that  are  commonly  used  for  the  strength  prediction 
of  composites,  it  can  be  concluded  that  fiber-based  failure 
criteria  (strength  of  materials  approach)  cannot  predict  size 
effects.  In  addition,  the  strength  of  materials  approach 
always  underpredicts  the  strength  of  notched  composites, 
with  errors  as  high  as  —49.1%  for  a  specimen  with  a 
2  mm  hole  diameter. 

The  Linear  Elastic  Fracture  Mechanics  approach  using 
a  hole  diameter  of  6  mm  for  calibration  predicts  the  size 
effect  accurately  for  specimens  with  hole  diameters  between 
6  mm  and  10  mm.  However,  Linear  Elastic  Eracture 
Mechanics  should  not  be  used  for  the  strength  prediction 
of  specimens  with  hole  diameters  equal  to  or  less  than 
2  mm,  or  for  larger  specimens  whose  failure  stresses  tend 
to  a  constant  value. 

The  point  stress  and  inherent  flaw  models  are  simple 
approaches  that  do  not  require  complex  EE  implementa¬ 
tions  yet  provide  reasonable  predictions  for  the  range  of 
hole  diameters  tested.  However,  the  accuracy  of  these  mod¬ 
els  relies  upon  the  measurement  of  the  characteristic  dis¬ 
tance  and  length  of  the  inherent  flaw  for  each  lay-up  and 
stacking  sequence. 

For  the  problems  selected,  the  continuum  damage 
model  proposed  predicts  with  good  accuracy  hole  size 
effects  in  composite  laminates  subjected  to  tension.  The 
model  requires  material  properties  that  are  measured  at 
the  ply  level  and  fracture  energies  that  are  measured  using 
both  standard  test  methods  and  novel  compact  tension  and 
compact  compression  test  methods.  The  continuum  dam¬ 
age  models  provides  not  only  the  final  failure  load,  but  also 
information  concerning  the  integrity  of  the  material  during 
the  load  history.  Furthermore,  the  finite  element-based 
damage  model  can  be  applied  to  structures  and  compo¬ 
nents  of  arbitrary  configurations  where  analytical  solutions 
could  not  be  developed. 
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Abstract 

In  this  article,  finite  element  analyses  of  mechanically  fastened  double-lap  joints  in  carbon/epoxy 
laminates  are  performed  using  a  progressive  damage  model  available  in  the  commercial  software 
ABAQUS.  An  alternative  damage  model,  implemented  into  a  VUMAT  user  subroutine,  is  also 
presented.  Two  failure  modes  are  considered:  catastrophic  net-section  tension-failure  and  non- 
catastrophic  accumulation  of  bearing  damage.  A  three-dimensional  mesh  is  used  for  the  analysis  and  in 
addition  to  results  for  static  implicit  analysis;  a  method  for  explicit  simulation  of  quasi-static  tests  is 
presented.  The  simulation  results  are  compared  with  experimental  data.  By  comparison  of  the  two 
damage  models  for  the  tension-failure  simulation,  it  can  be  shown  that  the  shape  of  the  damage  evolution 
law  for  fiber-tension  damage  is  perhaps  more  critical  than  the  fracture  energy  value.  Results  for 
simulation  of  bearing  damage  using  the  commercial  damage  model  are  presented  and  some  limitations  of 
the  model  are  discussed. 

1  Introduction 

With  the  increasing  use  of  fiber-reinforced  plastics  (FRPs)  in  aerospace  structures,  the 
analysis  of  mechanically  fastened  joints  in  composite  materials  has  become  a  key  aspect 
in  the  design  process.  It  is  well  known  that  mechanically  fastened  joints  perform  better 
in  metals  than  they  do  in  composite  structures.  The  joint  efficiency  in  a  metal  structure 
is  70%  -  80%  compared  to  40%  -  50%  in  a  composite  [1].  Some  reasons  for  the 
relatively  low  performance  of  bolted  composite  joints  are:  the  brittleness  of  the 
composite  material,  which  allows  little  stress  relief  around  the  loaded  hole;  material 
anisotropy,  leading  to  high  stress  concentrations;  and  low  through-thickness  strength  of 
classic  unidirectional  laminates,  causing  interlaminar  delamination.  Despite  these 
disadvantages,  mechanically  fastened  composite  joints  are  widely  used,  since  they 
provide  a  fast  and  efficient  way  of  substructure  assembly.  Due  to  the  complex  failure 
mechanisms,  their  design  however  relies  heavily  on  experiments  combined  with  semi- 
analytical  methods  [2].  If  it  is  possible  to  obtain  part  of  the  mechanical  properties 
needed  during  the  design  phase  via  numerical  analysis,  significant  cost  savings  can  be 
achieved.  Analysis  using  progressive  damage  models,  able  to  capture  the  physics  of  the 
failure  mechanisms  occurring  at  damage  initiation  and  damage  evolution  leading  to 
ultimate  failure  has  therefore  received  significant  attention  in  recent  years. 

In  general,  two-dimensional  finite  element  modelling  is  sufficient  for  the  majority  of 
linear  composite  laminate  analysis.  While  this  is  computational  efficient  and  preferable 
for  most  applications,  a  three-dimensional  model  may  be  suited  better  for  the  analysis  of 
a  bolted  composite  joint  in  a  quasi-isotropic  laminate.  In  a  3D-model,  cohesive  zone 
elements  can  be  included  to  capture  delamination  failure;  unsymmetrical  loading  of  the 
bolt  hole  (single-lap  joints)  can  be  considered  and  through-thickness  stresses  (clamping 
forces)  which  are  known  to  have  a  significant  effect  on  the  initiation  of  bearing  damage 
may  be  considered  [3]. 


2  Progressive  Damage  Models  for  Unidirectional  FRP  Laminates 

Two  progressive  damage  models  for  FRP  unidirectional  laminates  are  applied  in  the 
present  work.  The  first  model  recently  became  available  in  the  commercial  finite 
element  code  ABAQUS/Standard  6.6.1  and  ABAQUS/Explicit  6.7.1,  and  will  therefore 
be  referred  to  as  the  Abaqus-Model.  The  Hashin-criteria  is  used  for  damage  initiation  in 
this  model  [4], [5].  The  influence  of  damage  on  the  constitutive  material  model  is  based 
on  the  work  of  Matzenmiller  et  al.  [6]  and  damage  evolution  for  all  failure  modes  is 
governed  by  a  simple  linear  formulation,  used  by  Camanho  and  Davila  for  cohesive 
elements  [7].  A  detailed  description  of  the  Abaqus-Model,  including  its  numerical 
implementation,  is  presented  in  [8].  An  alternative  damage  model,  based  on  the  work  of 
Maimi  et  al.  [9]  is  also  used  in  this  study.  The  model  can  be  used  for  finite  element 
analysis  in  Abaqus/Explicit  via  a  VUMAT  user  subroutine  and  will  therefore  be 
referred  to  as  the  VUMAT-Model.  Maimi  applies  a  combination  of  the  EaRC03  and 
EaRC04  criteria  for  damage  initiation  [10], [11].  Rather  than  using  linear  softening, 
exponential  damage  evolution  laws  are  applied  to  describe  the  softening  response  for  all 
failure  modes  except  fiber  tension.  Eor  unidirectional  carbon/epoxy  laminates,  such  as 
the  material  used  in  this  study,  the  propagation  of  a  crack  perpendicular  to  the  fiber- 
direction  under  tensile  loading  can  be  divided  into  two  phases.  An  initial  and  rather 
brittle  fiber-matrix  failure  mechanism,  followed  by  a  tougher  fiber-bridging  and  fiber 
pull-out  phase  acting  at  a  lower  stress  level  [12].  To  account  for  the  different  damage 
mechanisms,  a  linear-exponential  law  is  therefore  used  for  the  fiber-tension  mode 
(Eigure  1,  b).  Eor  both  models,  the  area  under  the  stress-strain  curve  is  equal  to  the 
dissipated  fracture  energy  divided  by  a  characteristic  length  of  the  finite  element. 
References  [8]  and  [9]  provide  further  information  on  the  determination  of  the 
characteristic  element  length.  In  case  of  the  VUMAT-Model,  the  fracture  toughness 
determined  for  fiber-tensile  fracture  (Table  2)  is  divided  in  two  parts,  associated  with 
the  linear  and  exponential  softening  law.  In  addition  to  the  tensile  strength  ,  a  value 
representing  the  fiber  pull-out  strength  must  be  specified  for  the  VUMAT-model. 


(b)  VUMAT-model  [9] 


Eigure  1.  Damage  evolution  laws  for  fiber- tension 


3  Experiments 

Two  double-shear  bolted  joint  specimens  tested  in  the  context  of  developing  a  design 
methodology  for  mechanically  fastened  composite  joints  [2]  were  selected  in  this  study. 
The  specimen  geometry  and  dimensions  are  shown  in  Eigure  2  and  were  designed  to 
either  promote  pure  tensile-  or  bearing-failure.  Both  specimens  were  made  of  UD 
carbon/epoxy  prepreg  Hexed  IM7-8552  with  a  quasi-isotropic  lay-up  of 


[90/0/+45/-45]4s.  a  6  mm  steel  bolt  was  used  and  a  washer  with  an  outer-diameter  of 
12  mm  was  plaeed  on  either  side  of  the  laminate.  The  torque  applied  to  the  bolt 
eorresponds  to  a  finger-tight  assembly.  Surfaee  strain  was  measured  aeeording  to  the 
strain  gauge  positions  speeified  in  Figure  2.  Both  speeimens  were  tested  in  a 
eonventional  load  frame  at  a  quasi-static  displacement-rate  of  2  mm/min. 


Failure  Mode 

t 

1 

W 

d 

e 

ll 

l2 

Bearing 

3 

215 

36 

6 

18 

9 

50 

T  ension 

3 

200 

12 

6 

24 

15 

50 

h  h 


Figure  2.  Specimen  geometry  (dimensions  in  mm) 


4  FE  Model 

The  finite  element  model  used  for  both  specimen  types  will  be  explained  on  the  basis  of 
the  tension  failure  specimen  shown  in  Figure  3.  For  the  bearing  failure  specimen  a 
similar  mesh  was  used.  Due  to  the  lay-up  symmetry,  only  half  of  the  laminate  was 
modelled  and  symmetry  boundary  conditions  were  applied  at  the  laminate  symmetry 
plane.  One  element  per  ply  was  used  for  the  laminate  mesh,  which  was  divided  into  a 
coarse  mesh  area  away  from  the  hole  and  a  refined  mesh  area  around  the  hole  and  in  the 
direction  of  loading.  Both  mesh  regions  are  connected  via  a  TIE  constraint  which  is  a 
convenient  way  for  mesh  transition  as  opposed  to  a  paved  mesh  or  multi-point 
constraint  (MFC).  In  case  of  the  bolt,  only  the  length  of  the  shaft  in  contact  with  the 
hole  was  modelled.  Similar  symmetry  boundary  conditions  were  applied  to  the  nodes 
lying  in  the  laminate  symmetry  plane.  The  washer  is  accounted  for  by  a  distributed  load 
corresponding  to  the  bolt  torque  and  applied  to  a  surface  approximately  equal  to  the 
surface  area  of  the  washer.  Strain  in  the  loading  direction  is  obtained  from  two  element 
sets  in  the  first  layer  of  elements,  representing  the  90°  outer-ply,  at  the  strain  gauge 
position  of  the  test  specimen  (Figure  2).  Although  not  applied  in  the  present 
simulations,  the  3D  finite  element  mesh  was  developed  for  the  use  of  cohesive  zone 
delamination  elements  and  a  full  3D-formulation  of  the  VUMAT  damage  model,  which 
is  yet  to  be  implemented  into  an  Abaqus  subroutine.  For  the  Abaqus-Model,  where  two 
formulations  of  the  Hashin-criteria  are  available,  the  formulation  proposed  by  Hashin 
and  Rotem  was  selected  [4].  The  in-situ  effect  was  considered  for  both  damage  models. 
It  is  characterised  by  higher  transverse  tensile  and  shear  strengths  for  a  ply  constrained 
by  plies  with  different  fiber  orientations,  compared  to  the  strengths  of  the  same  ply  in  a 
unidirectional  laminate  [13].  For  the  tension- failure  specimen,  simulations  were 
conducted  using  the  implicit  and  explicit  Abaqus-model  as  well  as  the  VUMAT-model. 
The  bearing  failure  specimen  was  simulated  using  the  implicit  and  explicit  Abaqus- 
model.  Depending  on  the  damage  model  and  solver,  different  elements  were  used  for 
the  different  regions  of  the  finite  element  model.  The  selected  elements  are  summarised 
in  Table  1  where  SC8R  stands  for  a  reduced  integration  continuum  shell  element, 
similar  to  a  standard  solid  but  with  a  kinematic  and  a  constitutive  behaviour  similar  to  a 
conventional  shell.  The  Abaqus-Model  is  limited  to  elements  with  plane-stress 
formulation,  therefore  only  the  SC8R  element  can  be  used  in  a  3D-mesh.  C3D8  and 
C3D8R  represent  standard  solid  elements  in  a  fully  integrated  or  reduced  integration 


formulation,  respectively.  For  reduced  integration  elements,  default  hourglass  control 
parameters  were  selected.  In  case  of  the  implicit  Abaqus-model,  viscous  regularisation 
(VR)  had  to  be  used  to  obtain  a  converging  solution.  The  VR-parameters  were  selected 
according  to  a  similar  example  given  in  [8]. 

Table  1.  Finite  element  selection 


Solver 

Type 

Version 

Damage 

Model 

Laminate  Fine 

Mesh  Area 
Laminate  Coarse 

Bolt 

Abaqus/Standard 

6.6.1 

Abaqus 

SC8R 

C3D8 

C3D8 

Abaqus/Explieit 

6.7.1 

Abaqus 

SC8R 

C3D8R 

C3D8R 

Abaqus/Explieit 

6.6.1 

VUMAT 

C3D8R 

C3D8R 

C3D8R 

The  joint  is  loaded  via  a  velocity  boundary  condition  applied  to  a  selected  node-set  of 
the  bolt  mesh.  In  case  of  the  implicit  Abaqus-model,  this  velocity  corresponds  to  the 
actual  test  speed.  To  obtain  a  simulation  time  suited  for  an  explicit  simulation,  two 
modifications  were  applied  to  the  explicit  model.  The  test  speed  was  increased  by  a 
factor  of  1000  and  the  mass  density  was  scaled  by  a  factor  of  100,  resulting  in  a  10-fold 
increase  of  the  stable  time  increment. 


TIE  constraint 


laminate 
coarse  mesh 


washer  surface 


bolt 


strain  gauge  1 


laminate 
fine  mesh 


laminate  symmetry  plane 


Figure  3.  Finite  element  model  of  tension  failure  specimen 

To  evaluate  if  the  modified  FE-model  produces  an  appropriate  quasi-static  response,  the 
energy  balance  was  studied  [14]; 

El  +  Ey  +  E Pin  +  Epp  -  Epy  =  total  ~  const.  (4.1) 

El  is  the  internal  energy,  Ey  is  the  energy  absorbed  by  viscous  dissipation,  Ey^  is  the 
energy  absorbed  by  frictional  dissipation,  E^p  is  the  kinetic  energy,  E^  is  the  work  of 
external  forces  and  Ep^p^p  is  the  total  energy  in  the  system.  For  a  quasi-static  analysis, 
Efy  should  be  approximately  equal  to  Ei ,  while  as  Ey ,  Epp ,  E^p  and  Ep^p^p  should  be 

near  zero.  According  to  [14]  however,  a  5%  to  10%  value  of  the  kinetic  energy 
compared  to  the  internal  energy  is  acceptable.  A  graphical  representation  of  the  energy- 
balance  is  illustrated  in  Figure  4.  With  the  exception  of  a  slight  increase  in  the  total 
energy  Ep^p^p,  the  above  conditions  for  a  quasi-static  analysis  hold  for  the  entire 
duration  of  the  simulation.  The  total  energy  increase  prior  to  ultimate  failure  is  caused 
by  an  increase  of  Ey  and  .  These  two  energies  are  introduced  to  stabilise  the 


element  during  damage  evolution.  As  they  remain  small  eompared  to  and  ,  it 

was  coneluded  that  the  modifications  to  the  explicit  model  are  valid  and  the  simulation 
is  in  fact  representing  a  quasi-static  test. 


Figure  4.  Energy  balance,  tension-failure  simulation  (VUM AT -model) 

The  material  properties  used  for  both  progressive  damage  models  are  summarised  in 
Table  2.  Elastic  and  unidirectional  ply  strengths  can  be  obtained  from  standard  test 
methods  and  the  in-situ  ply  strengths  can  be  calculated  according  to  [13].  For  the 
fracture  energies,  a  test  standard  exists  only  for  matrix  tension  [15].  The  values  for 
tensile  and  compressive  fiber  fracture,  can  be  obtained  from  compact  tension  (CT)  and 
compact  compression  (CC)  tests  as  proposed  by  Pinho  et  al.  [16].  Matrix  compression 
fracture  energy  can  be  obtained  from  mode  II  end-notched  flexure  tests  (ENF)  and  a 
formulation  for  specified  in  [9].  According  to  this  formulation,  the  value  depends 
on  the  laminate  stacking  configuration.  In  the  present  study  the  value  for  a  strongly 
confined  laminate  is  used.  The  additional  parameters  used  for  the  VUMAT-model  are 
associated  with  the  different  damage  evolution  law  (compare  Figure  1,  b)  and  since  the 
LaRC  damage  initiation  criteria  considers  the  fracture  angle  for  compressive  transverse 
load  cases,  a  representative  fracture  angle  a  must  also  be  given. 

Table  2.  IM7-8552  material  properties 


(a)  Elastic  ply  properties  (b)  UD  ply  strength  (c)  In-situ  ply  strength  [MPa] 


El  171.42 

GPa 

2226.2  MPa 

Ply  configuration 

IS 

^is 

E2  =  E^  9.08 

GPa 

1200.1  MPa 

thin  outer 

101.4 

107.0 

Gi2  =  Gi3  5.29 

GPa 

62.3  MPa 

thin  embedded 

160.2 

130.2 

G23  3.98 

GPa 

yC 

199.8  MPa 

thin  embedded  (2t) 

113.3 

107.0 

*^12  “  *^13  0.32 

- 

SL 

92.3  MPa 

V23  0.5 

- 

(d)  Fracture  energies  [kJ/m^] 

(e)  Additional  VXJMAT-model  properties 

fiber  tension 

Gi, 

81.5 

fiber  pull-out  strength  X pQ 

232.3 

MPa 

fiber  compression 

Gi- 

106.5 

fracture  angle 

a 

53 

0 

matrix  tension 

G2. 

0.2774 

Gi^  ,  linear  softening  g/; 

31.5 

kJ/m^ 

matrix  compression 

G2- 

5.62 

Gi^  ,  exponential  softening 

50.0 

kJ/m^ 

5  Simulation  Results 


5 . 1  Tension-failure  simulation 

The  maximum  load,  ,  obtained  from  the  simulations  and  experiment,  is  summarised 

in  Table  3  and  Figure  5  shows  the  load-strain  response  at  strain  gauge  position  1  and  2. 
It  can  be  seen  that  both  simulations  using  the  Abaqus-model  significantly  overestimate 
the  ultimate  load  while  as  the  results  of  the  VUMAT-model  correlate  well  with  the 
mean  average  load  maximum  obtained  from  the  experiment.  Another  difference 
between  the  simulation  results  can  be  noticed  when  plotting  the  fiber  damage  parameter 
for  a  0°-ply  (Figure  6).  For  both  Abaqus-model  simulations,  the  crack  propagates  at  an 
angle  of  about  45°  to  the  fiber-direction  and  hence  follows  the  matrix  damage 
developing  in  the  neighbouring  45°-ply.  The  crack  in  the  VUMAT-model  localises  in  a 
plane  perpendicular  to  the  fibers,  as  observed  in  the  experiment. 

It  should  be  noted  that  various  damping  mechanisms,  such  as  mass-  or  stiffness- 
proportional  Rayleigh-damping  and  bulk  viscosity,  exist  for  the  Abaqus-model  while  as 
these  mechanisms  did  not  have  a  major  effect  on  the  user  material  specified  in  the 
VUMAT  subroutine.  Therefore  oscillations  measured  at  strain  gauge  1  could  not  be 
avoided  in  this  case. 


Table  3.  Maximum  load,  tension  failure  specimen 


^max  IkN] 

Experiment  (mean  average) 

9.477 

(minimum) 

9.232 

(-2.6%) 

(maximum) 

10.135 

(-r6.9%) 

Abaqus  Model,  implicit 

12.833 

(-r35.4%) 

Abaqus  Model,  explicit 

13.122 

(-r38.5%) 

VUMAT  Model 

9.454 

(-0.2%) 

Figure  5.  Load-strain  response,  tension- failure  specimen 


Note;  For  the  Abaqus-model,  a  separate  damage 
parameter  exists  for  fiber-tension  and  fiber- 
eompression  damage.  In  ease  of  the  VUMAT- 
model,  fiber-tension  and  fiber-eompression 
damage  are  combined  in  one  parameter. 
Therefore,  diagram  (c)  also  shows  fiber- 
compression  damage  at  the  bolt-hole  contact 
interface. 


(c)  VUMAT-model 


Figure  6.  Fiber  damage  in  layer  2  (0°)  at  maximum  load,  tension  failure  speeimen 


5.2  Bearing-failure  simulation 

Other  than  in  the  ease  of  the  tension-failure  speeimen,  where  ultimate  failure  is  elearly 
defined  by  the  maximum  load,  bearing-failure  is  a  non-eatastrophie  damage  mode, 
eharaeterised  by  a  progressive  accumulation  of  damage  and  permanent  hole  deformation 
[2],  As  a  result,  different  definitions  may  be  used  for  defining  bearing  strength  such  as 
the  onset  of  nonlinearity  or  the  bearing  strength  at  2%  bearing  strain  offset.  Figure  7  (a) 
shows  the  load-strain  response  obtained  from  simulation  and  a  representative 
experiment  at  strain  gauge  position  1 .  Diagram  (b)  illustrates  the  bearing  stress-bearing 
strain  curve  as  defined  by  the  ASTM  test  method  [17].  The  difference  between  the 
initial  bearing  stress  slope  of  simulation  and  experiment  can  be  explained  by  the 
different  method  of  obtaining  the  hole  elongation.  In  case  of  the  simulation,  the 
elongation  was  measured  directly  on  the  hole,  while  as  for  the  test  a  LVDT  was  attached 
to  the  test  rig  and  laminate,  similar  to  the  illustration  in  Fig.  10  (a)  of  [17]. 


(a)  Load-Strain  response  at  strain  gauge  1 


(b)  Bearing  stress  -  bearing  strain 


Figure  7.  Load-strain  and  bearing  stress-bearing  strain  response 


Table  4  summarises  the  bearing  strengths  defined  at  the  onset  of  nonlinearity  cr*',  and 

at  2%  bearing  strain  offset  for  the  implied  and  explieit  Abaqus-model  simulation 

and  the  experiment.  It  is  elear  from  the  experimental  results  that  the  2%  offset 
definition  is  assoeiated  with  signifieant  data  seatter  whieh  eomplieates  a  eomparison 
between  simulation  and  experiment.  It  was  shown  by  Camanho  and  Lambert  that 
damage  at  2%  strain  offset  has  progressed  to  a  state  of  through-the-thiekness  eracks 
spanning  several  plies  [2].  This  damage  state  ean  not  be  eaptured  by  the  present  finite 
element  model  and  therefore  the  onset  of  nonlinearity  (at  5%  deerease  of  the  initial 
ehord  modulus),  is  used  for  a  comparison  of  simulation  and  experiment  rather  than  the 
2%  offset  definition.  For  both  simulations,  the  predicted  bearing  strength  is  below  the 
experimental  value.  Figure  8  shows  the  extent  of  predicted  fiber-compression  damage 
at  2%  offset  bearing  strength. 


Table  4.  Bearing  strength 


—hr 

^onl 

[MPa] 

hr 

^2%off 

[MPa] 

Experiment  (mean  average) 

747 

870 

(minimum) 

738 

(-1.2%) 

747 

(-14.1%) 

(maximum) 

753 

(-rO.8%) 

958 

(-rlO.1%) 

Abaqus-model,  implieit 

645 

(-13.6%) 

851 

(-2.2%) 

Abaqus-model,  explieit 

610 

(-18.3%) 

689 

(-20.8%) 

(a)  Abaqus-model,  implicit  static 


(b)  Abaqus-model,  explicit 


Figure  8.  Fiber  damage  in  layer  2  (0°)  at  2%  offset  bearing  strength,  bearing  specimen 


6  Discussion  and  Conclusion 


6.1  Tension-failure  simulation 

In  Section  5.1  it  was  shown  that  the  ultimate  load  for  the  tension- failure  specimen  was 
significantly  over-predicted  by  the  Abaqus-model  (Table  3).  Since  the  same  material 
properties  were  used  for  all  simulations  and  viscous  regularisation  was  not  specified  for 
the  explicit  analysis,  it  is  assumed  that  the  difference  between  the  Abaqus-  and 
VUMAT-model  is  associated  with  the  damage  evolution  law.  It  was  further  noticed 
that  the  crack  for  the  Abaqus-model  develops  in  a  plane  inclined  at  an  angle  of  45°  to 
the  fibers.  In  an  attempt  to  create  a  damage  evolution  shape  similar  to  that  used  in  the 
VUMAT-model,  the  fracture  energy  for  fiber-tension  (compare  Table  2,  d)  was  reduced 
by  50%.  With  this  modification,  the  over-prediction  was  reduced  to  7.8%  for  the 
implicit  and  to  18.5%  for  the  explicit  formulation  of  the  Abaqus-model.  Comparing 
Figure  6  and  Figure  8,  it  can  be  seen  that  the  crack  has  shifted  towards  a  plane 


perpendicular  to  the  fiber-direction  with  the  improvement  most  pronounced  for  the 
implicit  model.  It  is  therefore  concluded  that  the  shape  of  the  fiber-tension  damage 
evolution  law  is  more  critical  than  the  actual  fracture  toughness  value  and  that  the 
formulation  chosen  in  the  VUM AT -model  is  able  to  represent  the  damage  mechanisms 
occurring  in  the  fiber-tension  damage  mode. 


(a)  Abaqus-model,  implicit  static  (b)  Abaqus-model,  explicit 

Figure  9.  Fiber-damage  in  layer  2  (0°)  at  maximum  load  for  Abaqus-model  with 

modified  fracture  toughness 

6.2  Bearing-failure  simulation 

In  Section  5.2  it  was  shown  that  conservative  results  can  be  obtained  for  the  bearing 
strength  at  onset  of  nonlinearity  (ONL),  using  either  the  implicit  or  explicit  formulation 
of  the  Abaqus-model  (Table  4).  With  damage  initiation  occurring  at  a  similar  stress 
level  of  about  450-460  MPa,  the  bearing  strength  results  for  the  ONL-formulation  of 
bearing  strength  were  relatively  close.  The  difference  in  the  nonlinear  region  of  the 
bearing  stress  curve,  comparing  the  implicit  and  explicit  solution  in  Figure  7  is  not 
entirely  understood  yet  and  may  partly  be  attributed  to  the  influence  of  viscous 
regularisation  on  the  evolution  of  damage.  As  for  the  tension- failure  specimen,  viscous 
regularisation  had  to  be  used  for  the  implicit  static  simulation  in  order  to  obtain  a 
converging  solution  but  was  not  used  in  case  of  the  explicit  simulations.  Further,  it  is 
possible  that  the  masses  considered  in  the  explicit  simulation,  and  hence  the 
modifications  as  described  in  Section  4,  may  have  an  influence  on  the  fiber- 
compression  damage  mode;  although  this  did  not  seem  to  be  critical  for  the  tension- 
failure  simulation. 

It  is  concluded  that  the  Abaqus-model,  using  a  simple  maximum  strength  criteria  and 
linear  damage  evolution  law  for  fiber-compression,  is  able  to  predict  a  lower  bound  for 
ONL  bearing  strength.  In  reality  however,  the  damage  mechanism  for  fiber- 
compression  is  more  complicated  with  fiber-kinking  occurring  in  the  0°  plies.  Further, 
the  model  can  not  account  for  the  stabilising-effect  of  through-thickness  stresses  [3], 
which  is  not  critical  if  a  relatively  low  level  of  clamping  pressure  is  used,  but  may  lead 
to  a  significant  under-prediction  at  higher  torque-levels. 

Acknowledgements 

The  financial  support  by  the  Air  Force  Office  of  Scientific  Research,  Air  Force  Material 
Command,  USAF,  under  grant  number  FA8655-06-1-3072  is  acknowledged.  The  U.S. 
Government  is  authorised  to  reproduce  and  distribute  reprints  for  Governmental 
purpose  notwithstanding  any  copyright  notation  thereon. 


References 


[1]  McCarthy  M.A.,  “BOJCAS:  Bolted  Joints  in  Composite  Aireraft  Struetures”,  Air 
and  Space  Europe,  2001;  3 

[2]  Camanho  P.P.,  Lambert  M.,  “A  design  methodology  for  meehanioally  fastened 
joints  in  laminated  eomposite  materials”,  Composites  Science  and  Technology, 
2006;  66:  3004-3020 

[3]  Park  HJ,  “Effeets  of  staeking  sequenee  and  elamping  foree  on  the  bearing 
strength  of  meehanioally  fastened  joints  in  eomposite  laminates”,  Composite 
Structures,  2001;  53:  213-221 

[4]  Hashin  Z.,  Rotem  A.,  “A  fatigue  failure  oriterion  for  fiber  reinforced  materials”. 
Journal  of  Composite  Materials,  1973;  7:  448-464 

[5]  Hashin  Z.,  “Failure  oriteria  for  unidireotional  fiber  oomposites”.  Journal  of 
Applied  Mechanics,  1980;  47:  329-334 

[6]  Matzenmiller  A.,  Lubliner  J.,  Taylor  R.L.,  „A  oonstitutive  model  for  anisotropio 
damage  in  fiber-oomposites“.  Mechanics  of  Materials,  1995;  20:  125-152 

[7]  Camanho  P.P.,  Davila  C.G.,  “Mixed-mode  deoohesion  finite  elements  for  the 
simulation  of  delamination  in  eomposite  materials”,  NASA-TM-2002-21 1737, 
2002 

[8]  Lapozyk  L,  Hurtado  J.A.,  “Progressive  damage  modelling  in  fiber-reinforoed 
materials”.  Composites  Part  A,  2007;  38:  2333-2341 

[9]  Maim!  P.,  Camanho  P.P.,  Mayugo  J.A.,  Davila  C.G.,  “A  eontinuum  damage 
model  for  eomposite  laminates:  Part  II  -  Computational  implementation  and 
validation”.  Mechanics  of  Materials,  2007;  39:  909-919 

[10]  Davila  C.G.,  Camanho  P.P.,  Rose  C.A.,  “Failure  eriteria  for  FRP  laminates”. 
Journal  of  Composites  Materials,  2005;  29:  323-345 

[11]  Pinho  S.T.,  Davila  C.G.,  lannueei  F.,  Robinson  P.,  “Failure  models  and  eriteria 
for  FRP  under  in-plane  or  three-dimensional  stress  states  ineluding  shear  non¬ 
linearity”,  NASA-TM-2005-2 13530,  2005 

[12]  Davila  C.G.,  Rose  C.A.,  “Superposition  of  eohesive  elements  to  account  for  R- 
curve  toughening  in  the  fracture  of  oomposites”,  Abaqus  Users’  Conference, 
2008 

[13]  Camanho  P.P.,  Davila  C.G.,  Pinho  S.T.,  lannueei  F.,  Robinson  P.,  “Prediotion  of 
in  situ  strengths  and  matrix  oraoking  in  oomposites  under  transverse  tension  and 
in-plane  shear”.  Composites  Part  A,  2006;  37:  165-176 

[14]  ABAQUS  6.6  “Getting  started  with  ABAQUS”,  Abaqus  Inc.,  2006 

[15]  “Standard  test  method  for  mode  I  interlaminar  fraoture  toughness  of 
unidirectional  fiber-reinforced  polymer  matrix  oomposites”  ASTM  D  5528-01, 
Amerioan  Sooiety  for  Testing  and  Materials  (ASTM),  West  Conshohooken,  PA, 
USA 

[16]  Pinho  S.T.,  Robinson  P.,  lannueei  F.,  “Fraoture  toughness  of  the  tensile  and 
oompressive  fibre  failure  modes  in  laminated  oomposites”.  Composites  Science 
and  Technology,  2006;  66:  2069-2079 

[17]  “Standard  test  method  for  bearing  response  of  polymer  matrix  eomposite 
laminates”  ASTM  D  5961/D  5961  M-05,  Amerioan  Sooiety  for  Testing  and 
Materials  (ASTM),  West  Conshohooken,  PA,  USA 


Appendix  D:  Paper  published  in  the  Journal  of  Com¬ 
posite  Materials 


56 


A  three-dimensional  damage  model  for  transversely 
isotropic  composite  laminates 

P.  Maimf^,  P.P.  Camanho'^’*,  J.A.  Mayugo^ 

^AMADE,  Universitat  de  Girona,  Campus  Montilivi  s/n,  Girona,  Spain 

^ DEMEGI,  Eaculdade  de  Engenharia,  Universidade  do  Porto,  Rua  Dr.  Roberto  Erias,  4^00-465, 

Porto,  Portugal 


Abstract 

This  paper  proposes  a  fully  three-dimensional  continuum  damage  model,  developed  at  the  sub-ply 
level,  to  predict  in  an  integrated  way  both  the  intralaminar  and  the  interlaminar  failure  mechanisms 
that  occur  in  laminated  fiber-reinforced  polymer  composites.  The  constitutive  model  is  based  on  the 
assumption  that  the  composite  material  is  transversely  isotropic,  and  accounts  for  the  effects  of  crack 
closure  under  load  reversal  cycles.  The  damage  model  is  implemented  in  an  implicit  finite  element 
code  taking  into  account  the  requirement  to  ensure  a  mesh-independent  computation  of  the  dissipated 
energy.  The  comparison  between  the  model  predictions  and  published  experimental  data  indicates 
that  the  model  can  accurately  predict  the  effects  of  transverse  matrix  cracks  on  the  residual  stiffness 
of  quasi-isotropic  laminates,  the  interaction  between  transverse  matrix  cracks  and  delamination,  and 
final  failure  of  the  laminate. 

Key  words:  Fracture,  Damage  Mechanics,  FEA. 


[Table  1  about  here. 


1  Introduction 


It  has  been  widely  recognized  that  one  of  the  most  signihcant  barriers  to  the  increased  use  of 
composite  materials  is  the  inability  to  predict  accurately  structural  failure  [1].  The  prediction 
of  structural  failure  in  laminate  composites  is  particularly  challenging  when  both  delamination 
and  intraply  failure  mechanisms,  such  as  matrix  cracking  or  fiber  failure,  contribute  to  the 
fracture  process. 

Delamination  is  normally  simulated  using  methods  based  on  Linear-Elastic  Fracture  Mechan¬ 
ics,  such  as  the  Virtual  Crack  Closure  Technique  [2],  or  using  cohesive  formulations  [3]-[9].  The 
onset  of  intralaminar  failure  mechanisms  is  normally  predicted  using  ply-based  failure  criteria 
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11  [10]- [15].  Generally,  failure  criteria  alone  are  unable  to  predict  the  collapse  of  composite  struc- 

12  tures.  To  predict  failure  initiation,  propagation  and  final  collapse  it  is  necessary  to  combine 

13  the  ply-based  failure  criteria  with  appropriate  damage  models.  When  the  laminate  is  uniformly 

14  stressed,  and  when  transverse  matrix  cracks  are  the  main  failure  mechanism,  it  is  possible  to 

15  use  analytical  or  semi-analytical  solutions  to  relate  the  applied  load  to  the  residual  stiffness  and 

16  strength  of  the  laminate  [16]- [21].  However,  the  simulation  of  the  propagation  of  intralaminar 

17  failure  mechanisms  in  composite  structures  with  complex  geometries  requires  models  based  on 

18  the  formalism  of  continuum  damage  mechanics  [22]- [30]. 

19  There  are  several  relevant  structural  applications  of  laminated  composites  where  both  delami- 

20  nation  and  ply  failure  mechanisms  are  relevant,  interacting,  energy  dissipation  mechanisms.  For 

21  example,  in  composites  subjected  to  low  velocity  impact,  in  skin-stiffener  terminations  or  in 

22  ply-scaled  notched  laminates.  The  approach  normally  used  to  simulate  both  delamination  and 

23  ply  failure  is  to  combine  cohesive  elements  that  simulate  delamination  with  continuum  damage 

24  models  that  simulate  ply  damage  [31]-[33].  Although  this  mesomechanical  approach  has  proved 

25  to  be  successful  for  some  structural  conhgurations  [31],  there  are  some  fundamental  problems 

26  that  prevent  the  general  application  of  this  methodology  that  uses  two  different  kinematic  repre- 

27  sentations  for  interlaminar  and  intralaminar  failure  mechanisms.  For  example,  in  the  numerical 

28  simulation  of  the  interaction  between  transverse  matrix  cracks  and  delamination  it  is  necessary 

29  to  capture  the  high  stresses  at  the  tip  of  the  transverse  crack.  Using  a  mesomechanical  model, 

30  it  is  not  possible  to  capture  this  interaction  because  the  elements  where  the  transverse  crack 

31  is  predicted  soften  without  being  able  to  accurately  capture  the  stress  held  at  the  interface. 

32  Furthermore,  even  if  strain-softening  constitutive  models  are  used  in  mesomechanical  models 

33  to  predict  transverse  matrix  cracking  in  multidirectional  laminates,  the  hnite  element  where 

34  transverse  cracking  is  predicted  does  not  unload  the  adjoining  elements  that  represent  the  same 

35  ply,  and  therefore  is  unable  to  represent  a  transverse  crack. 

36  To  accurately  predict  the  interaction  between  intralaminar  and  interlaminar  failure  mechanisms 

37  it  is  essential  to  have  a  good  kinematic  representation  of  the  different  failure  mechanisms. 

38  This  has  been  realized  by  Wisnom  and  co-workers  [34], [35]  in  their  simulations  of  fracture  of 

39  notched  and  unnotched  specimens,  using  cohesive  zone  models  to  represent  both  intralaminar 

40  and  interlaminar  failure  mechanisms.  The  cohesive  zone  models  provide  an  accurate  kinematic 

41  representation  that  enabled  the  successful  simulation  of  complex  phenomena  such  as  size  effects 

42  in  both  ply-scaled  and  sublaminate  scaled  composite  notched  specimens  [34],  and  the  fracture 

43  of  unnotched  scaled  quasi-isotropic  specimens  [35].  Although  the  use  of  cohesize  zone  models 

44  has  proved  to  be  an  accurate  technology  to  simulate  the  interaction  between  failure  mechanisms 

45  there  are  two  main  limitations  to  its  use.  The  hrst  limitation  is  the  need  to  know  in  advance 

46  the  planes  of  crack  propagation.  There  are  several  situations  where  the  orientation  of  the  crack 

47  plane  is  not  know  a  priori.  For  example,  when  a  ply  is  subjected  to  both  transverse  compression 

48  and  in-plane  shear  the  fracture  plane  depends  on  the  relation  between  these  two  components 

49  of  the  stress  tensor  [12].  The  second  limitation  to  the  use  of  cohesive  zone  models  is  the  need 

50  to  introduce  cohesive  hnite  elements  at  every  single  interface  where  a  crack  may  develop. 

51  To  overcome  some  of  the  difficulties  in  the  simulation  of  the  interaction  between  failure  mecha- 

52  nisms,  the  objective  of  this  work  is  to  formulate  a  fully  three-dimensional  damage  model  at  the 

53  sub-ply  level  that  is  able  to  represent  both  interlaminar  and  intralaminar  failure  mechanisms 

54  without  previous  knowledge  of  the  orientation  of  the  failure  planes.  The  sub-ply  level  constitu- 

55  five  model  and  the  corresponding  computational  implementation  are  described  in  the  following 

56  sections.  The  model  is  then  validated  by  comparing  the  numerical  predictions  with  published 
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experimental  results. 


58  2  Damage  model  for  a  transversely  isotropic  laminate 


59  2.1  Constitutive  tensor 


60  Consider  a  transversely  isotropic  material  and  a  vector  ei  =  {1,0,0}^,  parallel  to  the  fiber 

61  direction.  In  the  transversely  isotropic  plane,  orthogonal  to  ei,  there  is  a  set  of  orthonormal 

62  vectors  {e2,  03}  that  define  a  plane  where  the  shear  strain  723  is  zero.  In  order  to  correctly  detect 

63  crack  closure  under  load  reversal  cycles  a  set  of  ortonormal  vectors  {61,62,63}  are  defined. 

64  Taking  e  =  (en,  622,  eaS)  2ei2, 2ei3,  2623}^  as  the  components  of  the  strain  tensor  in  the  global 

65  coordinate  system,  the  relation  between  the  strain  tensor  in  the  local  and  global  coordinate 

66  system  is: 


"  =Te 

0 


(1) 


67  The  transformation  matrix  T  relates  the  strains  in  the  material  coordinate  system  to  the  strains 

68  in  the  coordinates  defined  by  the  vectors  {61,62,63},  e  =  {^n,  £22,  £33, 7i2,  Without  loss 

69  of  generality,  it  is  possible  to  assume  that  the  one  direction  of  the  global  coordinate  system 

70  coincides  with  61.  Therefore  the  transformation  matrix  can  be  written  as: 
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1  0 

0  cos^  6 

0  sin^  e 

0  0 

0  0 


0 

sin2  0 
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0  0  0 
0  0  cos  9  sin  9 

0  0  —  sin  6^  cos  9 


0  cos  6*  —  sin  6*  0 

0  sin  9  cos  9  0 


0  —  2sin6'cos6' 2sin6*cos6*  0  0  cos^ 6^  —  sin^ 0 


(2) 


71  The  angle  9  is  determined  by  enforcing  that  the  shear  strain  723  is  zero,  i.e.,  tan  {—29)  = 

72  Having  defined  the  coordinate  system  that  is  the  basis  for  the  derivation  of  the  the  constitutive 

73  model,  it  is  now  necessary  to  propose  a  suitable  form  for  the  specific  free  energy.  Assuming 

74  a  constant  density,  the  total  complementary  free  energy  is  given  as  Jy  'ipdV ,  where  ip  is  the 

75  complementary  free  energy  per  unit  volume.  The  proposed  definition  for  the  complementary 

76  free  energy  per  unit  volume  is: 
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+  [/5i1(Tii  +  /322  ((T22  +  (T33)]  am 


77  where  0:11  and  0:22  are  the  coefficients  of  thermal  expansion  in  the  longitudinal  and  transverse 

78  direction,  respectively,  jdn  and  (322  the  coefficients  of  hygroscopic  expansion  in  the  longitudinal 

79  and  transverse  direction,  respectively.  AT  and  AM  are  the  differences  in  temperature  and 

80  moisture  content  with  respect  to  the  corresponding  reference  values,  cfi,  d2  and  d^,  are  the 

81  damage  variables.  The  strain  tensor  is  calculated  as: 


dib  d'^'ib 

e  =  —  =  H  ;  (T  +  q;AT  +  pAM,  with  H  =  7- - —  (4) 

da  da  0  da 

82  The  compliance  tensor,  H,  relates  the  elastic  strains  with  the  elastic  stresses.  This  tensor 

83  depends  on  the  value  of  damage  variables.  It  is  assumed  that  damage  is  represented  by  a  set 

84  of  variables  that  affect  the  longitudinal,  the  transverse  and  the  shear  modulus.  The  damage 

85  variables  related  to  longitudinal  and  transverse  directions  change  when  the  normal  stresses 

86  switch  from  positive  to  negative  or  vice-versa.  The  damage  variable  di  represents  cracks  in 

87  planes  normal  to  the  fiber  direction,  whereas  the  damage  variables  d2  represents  cracks  in 

88  planes  parallel  to  the  hber  direction.  The  damage  variable  d^  affect  the  shear  moduli  G12  and 

89  G13.  The  transverse  damage  variable  is  not  able  to  detect,  at  constitutive  level,  the  directionality 

90  of  cracks.  However,  the  directionality  of  cracks  is  detected  at  macroscopic  level  as  the  locus  of 

91  the  damaged  points.  From  (4),  the  compliance  tensor  is  dehned  as: 

1  ffi2 

—  di)  El  El 

1 

El  (1  —  d2)  E2 

_  ffi2  _  ^23 

El  E2  ( 1 

0  0 

0  0 

92  where  Ei  and  E2  are  the  longitudinal  and  transverse  Young  modulus,  respectively,  z/12  is  the 

93  major  Poisson  ratio  and  z/23  is  the  Poisson  ratio  in  the  transverse  isotropic  plane.  G12  is  the 

94  shear  modulus,  di,  d2  and  d^  are  the  damage  variables  in  the  directions  dehned  by  the  vectors 

95  (61,62,63).  These  damage  variables  depend  on  the  longitudinal  {dL±)  and  transverse  ((ir±) 

96  damage  variables  as: 
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97  where  (x)  is  the  McCauley  operator  dehned  as  {x)  ■.=  {x  +  \x\)  /2.  The  shear  damage  variable 

98  is  not  influenced  by  the  sign  of  the  shear  stress  components,  i.e.,  d^  =  ds-  The  shear  damage 

99  variable  is  influenced  by  the  normal  stresses  that  produce  friction  between  the  crack  faces 

100  allowing  stress  transfer  and  dissipation  under  shear  loads  [36].  This  effect  is  neglected  in  the 

101  present  model. 

102  It  should  be  noted  that  the  closure  effect  in  the  transversely  isotropic  plane  is  activated  inde- 

103  pendently  in  the  the  directions  62  and  63.  Therefore  the  material  becomes  orthotropic  when 

104  (T22  has  a  different  sign  than  that  of  (T33. 


105  The  coaxially  of  stresses  and  strains  in  the  transverse  isotropic  plane  is  enforced  and  the  cor- 

106  respondent  shear  modulus  is  evaluated  as:  G23  =  If  the  stresses  in  the  transversely 

107  isotropic  plane  have  the  same  sign,  and  if  the  damage  variables  have  the  same  value  (^2  =  d^), 

108  the  shear  modulus  is  given  as:  G23  =  2[i+t23(i-L)] '  important  to  note  that,  unlike  the 

109  majority  of  the  rotating  crack  models,  this  model  does  not  exhibit  a  negative  shear  modulus, 
no  which  is  a  physically  inadmissible  result  [37].  This  is  due  to  the  fact  that  the  damage  is  isotropic 

111  if  the  transverse  stresses  have  the  same  sign. 

112  The  stress  tensor  in  the  global  coordinate  system  (^)  is  calculated  as: 


113  2.2  Dissipation 

114  The  thermodynamic  forces  Ym  {M  =  1,  2,  3,  6)  are  calculated  as: 
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115  The  rate  of  dissipation  is  expressed  in  terms  of  the  thermodynamic  forces  and  damage  variables 

116  as: 
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117  Due  the  particular  form  of  the  complementary  free  energy  selected,  it  is  observed  that  the 

118  thermodynamic  forces  (Dm)  are  always  positive.  Therefore  the  condition  of  positive  evolution 

119  of  damage  variables  (c^m  >  0)  ensures  a  positive  energy  dissipation.  The  crack  closure  effect 

120  under  load  reversal  cycles  does  not  result  in  spurious  energy  dissipation  because  the  conjugated 

121  thermodynamic  forces  are  zero  when  such  a  loading  scenario  takes  place  [38]. 

122  2.3  Damage  activation  functions 

123  The  damage  activation  functions  dehne  the  elastic  domain  under  general  stress  states.  The 

124  elastic  domain  is  dehned  here  by  three  damage  activation  functions,  that  are  represented  by 

125  three  surfaces  in  the  strain  space. 

126  The  selection  of  the  damage  activation  function  depends  upon  the  different  failure  mechanisms 

127  of  the  material  system.  The  main  assumption  of  the  present  selection  of  damage  activation 

128  functions  is  that  the  shear  stresses  in  the  transversely  isotropic  plane,  (Ji2  and  ais,  create 

129  cracks  orientated  in  a  plane  with  a  normal  vector  in  the  (02, 63)  plane.  This  response  is  typical 

130  of  nnidirectional  composites  in  which  the  hbers  enforce  the  matrix  cracks  to  growth  along  their 

131  direction. 

132  The  damage  activation  fnnctions  Tjv  {N  =  L+,L— ,T)  are  dehned  as: 

Fl+  =  (t)L+  -  rL+  <  0 

Fl-  =  0L-  -  <  0  (10) 

Fj'  =  (pj'  —  Tj'  ^0 

133  where  Fl+  dehnes  the  elastic  domain  for  longitndinal  tensile  failure,  Fl-  dehnes  the  elastic 

134  domain  for  longitndinal  compressive  failnre,  and  Ft  dehnes  the  elastic  domain  for  transverse 

135  failure. 

136  The  loading  functions  {N  =  L+,  L—,T)  depend  on  the  strain  tensor,  and  on  the  elastic 

137  and  strength  properties.  The  internal  variables  r^v  {N  =  L+,  L—,T)  take  an  initial  value  of  1 

138  when  the  material  is  undamaged,  and  they  increase  with  damage.  The  internal  variables  of  the 

139  constitntive  model  are  related  to  the  damage  variables  dM  {M  =  L+,  L—,T+,T—,  S)  by  the 

140  damage  evolution  laws. 

141  2.4  Loading  functions 

142  A  simple  non  interacting  maximnm  strain  or  stress  criteria  resnlts  in  accurate  predictions  of 

143  the  onset  of  longitudinal  damage  of  polymer-based  composite  materials  nnder  tensile  stresses 

144  [14], [15].  The  maximum  strain  criterion  is  used  for  longitudinal  tensile  loading: 
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(11) 


4>l+  -  (^ii) 

145  Longitudinal  failure  of  unidirectional  composite  materials  under  compressive  stresses  is  a  far 

146  more  complex  phenomenon.  The  compressive  failure  is  the  result  of  a  complex  sequence  of 

147  damage  mechanisms  that  culminate  in  the  formation  of  a  kink  band. 

148  The  LaRC03-04  [14], [15]  failure  criteria  postulates  that  hber  kinking  is  triggered  by  onset  of 

149  damage  in  the  matrix.  Under  this  circumstance,  the  fibers  loose  their  lateral  support  and 

150  collapse  under  the  compressive  load.  The  failure  load  depends  on  the  initial  hbre  misalignment, 

151  and  on  the  rotation  of  the  fibers  as  a  function  of  the  applied  stresses.  In  this  model  a  simple 

152  non-interacting  maximum  strain  is  used: 

fe-  =  (12) 

153  It  is  clear  that  such  a  simple  criterion  is  unable  to  account  for  the  effects  of  the  shear  stresses 

154  on  hber  kinking.  The  development  of  a  more  accurate  failure  criterion  for  three-dimensional 

155  stress  states  will  be  addressed  in  future  work. 

156  The  transverse  loading  function  dehnes  the  onset  of  transverse  failure  mechanisms.  The  loading 

157  function  has  to  match  the  three  uniaxial  loads  that  produce  matrix  cracking:  transverse  tension, 

158  transverse  compression  and  longitudinal-transverse  shear.  The  proposed  loading  function  is: 


159  where  a  are  the  effective  stresses,  calculated  using  the  undamaged  stiffness  tensor,  given 

160  by  (5)  with  dj  =  0  (f  =  1,  2,  3,  6),  as  d  =  Hg  ^  :  £.  The  transverse  loading  function  has  the  same 

161  form  as  the  criteria  presented  by  Christensen  [39]. 

162  Figure  1  shows  the  transverse  damage  activation  function.  As  previously  mentioned,  the  inter- 

163  hber  shear  stress  produce  transverse  cracks.  Therefore,  the  transverse  damage  must  be  activated 

164  under  uniaxial  shear  loads.  Furthermore,  experimental  results  and  failure  criteria  developed 

165  following  the  Mohr-Coulomb  theory,  as  such  as  the  Puck  [12]  and  the  LaRC  [14]- [15]  criteria, 

166  demonstrate  that  moderate  values  of  transverse  compression  increase  the  shear  strength.  As 

167  shown  in  Figure  1  c),  this  ehect  is  accounted  for  in  the  proposed  loading  function. 

168  [Fig.  1  about  here.] 

169  2.5  Internal  variables 

170  Neglecting  viscous  ehects,  the  damage  activation  functions  must  be  negative.  If  Fat  <  0  the 

171  material  is  in  the  elastic  regime.  When  the  damage  activation  criteria  is  satished.  Fat  =  0,  it  is 

172  necessary  calculate  the  gradient  ^at-  If  07v  is  negative  or  zero,  the  state  is  one  of  unloading  or 
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173  neutral  loading,  respectively.  The  different  states  of  the  material  response  are  mathematically 

174  represented  by  the  Kuhn- Tucker  conditions,  >  0;  Fjq  <  0;  rjsqFpq  =  0.  If  07v  is  positive, 

175  damage  increases,  and  the  consistency  condition  has  to  be  satished,  i.e.,  Fm  =  0  ^  Fjq  =  0. 

176  If  the  internal  variables  are  exclusively  dependent  on  the  damage  variables,  and  if  the  loading 

177  functions  depends  on  the  strain  tensor,  the  constitutive  model  can  be  explicitly  integrated  [40] , 

178  [41],  Applying  the  consistency  condition,  the  internal  variable  tt  is  calculated  as: 

=  max  |l,  max  (14) 

179  The  evolution  of  the  longitudinal  elastic  domain  thresholds  for  tensile  or  compressive  stresses 

180  are  coupled.  The  elastic  domain  threshold  defines  the  level  of  elastic  strains  that  can  be  attained 

181  before  the  accumulation  of  additional  damage. 

182  Under  longitudinal  tensile  stresses,  the  fracture  plane  is  generally  perpendicular  to  the  hber 

183  direction.  When  reversing  the  load,  the  cracks  close  and  can  still  transfer  load.  However,  the 

184  broken  and  misaligned  hbers  do  not  carry  any  additional  load.  Therefore,  the  compressive 

185  stiffness  is  influenced  by  longitudinal  damage.  However,  the  elastic  domain  is  assumed  to  remain 

186  unchanged.  Under  longitudinal  compression,  damaged  material  consisting  of  broken  hbers  and 

187  matrix  cracks  forms  a  kink  band,  and  there  is  not  a  unique  orientation  for  the  damage  planes. 

188  When  the  loads  are  reversed,  the  cracks  generated  in  compression  open  and  the  elastic  domain 

189  threshold  increases.  To  represent  these  phenomena,  the  evolution  of  the  longitudinal  internal 

190  variables  is  dehned  as: 

rL+  =  0L+  and  tl-  =  0 

;  ,  .  0L-  if  rL+  <  Tl- 

tl-  =  4>l-  and  rL+  =  < 

[  0  if  rL+  >  tl- 

191  The  integration  of  the  previous  expressions  results  in: 

=  max  |l,  aiax 
rL_  =  max  |l,  max 

192  2.6  Damage  evolution  laws 

193  The  dehnition  of  the  damage  evolution  laws,  which  relate  the  internal  variables  with  the  damage 

194  variables  dMir^),  is  required  to  fully  dehne  the  constitutive  model. 

195  When  the  material  is  undamaged,  the  internal  variables  r^r  take  the  initial  value  of  1,  and 

196  dMi^N  =  1)  =  0.  Equations  (14)  and  (15)  dehne  the  evolution  of  the  internal  variables  ensuring 
that  tn  >  0.  As  shown  in  equations  (8)  and  (9),  the  condition  for  positive  dissipation  is 
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Tension  loading: 
Compression  loading: 
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198  satisfied  if  cIm  >  0.  The  condition  for  positive  dissipation  is  automatically  fulfilled  if  the  damage 

199  evolution  law  satishes  the  condition  ddM/dr^  >  0.  When  the  material  is  completely  damaged, 

200  a  fracture  plane  is  created,  the  strains  are  localized  in  a  plane  in  which  r^r  — >  cxd  and  the  related 

201  components  of  the  stiffness  tensor  are  zero,  dMifN  — >  cio)  =  1. 

202  The  evolution  of  internal  variables  can  result  in  material  hardening  or  softening  depending  the 

203  damage  law.  If  the  stress-strain  response  result  in  a  softening  relation,  the  deformations  localize 

204  in  a  plane,  and  a  localization  limiter  has  to  be  introduced  in  the  model  to  correctly  compute 

205  the  energy  dissipated. 

206  The  procedure  followed  in  this  model  to  ensure  a  correct  computation  of  the  energy  dissipated 

207  is  based  on  the  Crack  Band  Model  proposed  by  Bazant  [42].  Using  equation  (9)  it  is  possible 

208  calculate  the  dissipated  energy  under  an  uniaxial  test  as: 

9m  =  rYMdMdt=  r  Yu^dr^  =  M  =  L+,L-,T+,S  (16) 

JO  dr'N  t 

209  where  i  is  the  characteristic  length  of  the  hnite  element,  and  Gm  is  the  associated  fracture 

210  toughness.  The  fracture  toughness  Gl+  corresponds  to  a  crack  that  propagates  in  a  plane 

211  perpendicular  to  the  hber  direction  under  mode  I  loading  and  the  fracture  toughness  Gl-  is 

212  related  to  hbre  kinking.  There  are  no  standard  test  methods  to  measure  these  properties.  How- 

213  ever,  Pinho  and  co-authors  [43]  developed  new  compact-tension  (CT)  and  compact-compression 

214  (CC)  test  methods  that  can  be  used  to  measure  Gl+  and  Gl_.  The  fracture  toughness  Gt+ 

215  and  Gs  correspond  respectively  to  matrix  cracking  for  mode  I  and  mode  II  loading.  These 

216  properties  can  be  measured  using  standard  double-cantilever  beam  and  end-notched  flexure 

217  test  specimens. 

218  If  the  characteristic  element  size  {£)  is  greater  than  a  critical  value  the  material  response  results 

219  in  snap-back,  and  the  energy  dissipation  would  be  overpredicted.  To  prevent  this  problem,  the 

220  characteristic  element  size  must  be  lower  that  a  critical  value  given  by: 

i  =  ^  ;  M  =  L+,L-,T+,S  (17) 

221  If  the  element  is  larger  that  the  maximum  size  prescribed,  and  a  mesh  rehnement  is  unfeasible, 

222  the  snap-back  in  the  constitutive  model  can  be  avoided  by  reducing  the  corresponding  strength 

223  according  to  [44]: 

=  M  =  L+,L-,T+,^  (18) 

224  2.6.1  Transverse  tension  damage  law 

225  Under  transverse  tension,  damage  localizes  in  a  fracture  plane  without  any  previous  inelastic 

226  material  behavior.  The  linear  softening,  or  cohesive,  law  shown  in  Figure  2a),  is  proposed. 

227  [Fig.  2  about  here.] 
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228 


For  a  linear  softening  law,  the  stress  strain  response  is  given  as:  (T22  =  HT+E2S22  +  Y  01 
229  (T22  =  (1  ~  dT+)E2e22-  Using  the  resulting  £22  in  (13)  and  (14),  the  damage  law  is  calculated  as: 


dT+  —  1  —  — 


‘lAj'-y-  (1  —  Yj' 


(19) 


230  where  E[t+E2  is  the  incremental  stiffness  under  uniaxial  stress,  and  the  remaining  parameters 

231  are: 


At+  — 


(1  +  ^2zHT+f‘ 

(1  +  ^^23)^ 


(Yc  —  Yt)  (1  —  EIt+  (2z/i2Z2’21  +  ^’23))  ,  o  (^  +  ^23Ht+)  1^23  (1  ”  Ht+)  Yt 
Et+  = - ; - ^ - h  2- 


1  -  Z/23  -  2z/i2Z/21 


(1  +  ^^23)" 


Ct+  — 


1223  (1  —  Ht+)  Yt\  (Yc  —  Yt)  (1  —  Ht+)  Yt  {1223  +  ‘^1212^21) 


1  +  ^^23 


1  —  ^23  ~  2z/i2Z22i 


232  The  parameter  Ht+  is  calculated  by  applying  the  crack  band  model,  equation  (16): 


Ht+  = 


Y^l 


Y:^t 


2GT+E2 


<  0 


(20) 


233  2.6.2  Transverse  compression  damage  law 

234  The  law  proposed  to  simulate  damage  evolution  under  transverse  compression  is  shown  in 

235  Figure  2  a),  and  it  is  given  as: 


Bt-  +  ijBl_  -  4:At-  (Ct-  -  YTYcrV) 

236  where  HT-E2  is  the  incremental  stiffness  under  an  uniaxial  compression  load.  Ht-  dehnes  the 

237  hardening  in  the  compressive.  The  parameters  At-,  Bt-  and  Ct-  are  determined  following 

238  the  procedure  outlined  in  the  previous  section  as: 


At— 

Bt— 

Ct— 


(1  +  ^23Ht-Y 
(1  +  ^^23)^ 

(Yc  —  Yt)  (1  —  Ht-  (212121221  +  1223))  _  2(^  +  ^‘23Ht-)  1^23  (1  —  Ht-)  Yc 
1  ~  ^23  ~  2z/i2Z/21  (1  +  1223)“^ 

1 1223  (1  -HT-)Ycy  ^  (Yc  -  Yt)  (1  -  Ht-)  Yc  (z/23  +  212,21^21) 

V  ^  +  ^23  )  1  -  z/23  -  2z/i2Z/21 
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239  2.6.3  Longitudinal  tension  damage  law 

240  Crack  propagation  in  the  longitudinal  direction  involves  different  energy  dissipating  mecha- 

241  nisms  such  as  fiber  fracture,  fiber-matrix  pull-out  and  matrix  fracture.  When  different  physical 

242  mechanisms  are  involved  in  crack  propagation,  a  two-part  damage  evolution  law  should  be  used 

243  [45].  To  accurately  represent  the  propagation  of  longitudinal  failure  mechanisms,  the  authors 

244  have  previously  proposed  a  softening  law  defined  by  the  linear-exponential  relation  shown  in 

245  Figure  2b)  [26],  [27].  It  is  observed  in  Figure  2b)  that  the  softening  response  is  linear  until  the 

246  stress  reaches  the  pull-out  stress,  Xpo,  and  the  corresponding  energy  dissipation  per  unit  area 

247  is  Gp_|_.  These  material  parameters  can  be  measured  using  a  recently  proposed  analysis  method 

248  for  the  resistance  curve  measured  in  the  compact  tension  test  specimens  [28]. 

249  As  the  strains  continue  to  increase,  the  softening  response  follows  an  exponential  law  and  the 

250  energy  dissipated  per  unit  area  is  Gf+.  The  definition  of  the  damage  evolution  law  used  here 

251  is  explained  in  detail  in  previous  papers  [26], [27],  and  the  resulting  equations  are: 


dL+  —  (1  +  Hi 


rL+ 


dL+  =  1 


Xi 


PO 


XTrL+ 


rL+ 


exp 


A 


rL+  -  rL+ 


if  rp+  <  rl^ 


if  rp+  > 


(22) 


252  where  the  parameters  Hl,  and  are  defined  in  [26], [27]. 

253  2.6.4  Longitudinal  compression  damage  law 

254  The  compressive  stiffness  is  influenced  by  both  the  damage  produced  under  compression  and 

255  under  tension.  If  the  material  is  damaged  in  tension  and  the  loads  are  reversed  until  the  material 

256  is  loaded  in  compression,  the  cracks  close  allowing  the  stress  transfer.  Although  the  cracks 

257  are  closed,  the  fibers  lose  the  alignment  and  do  not  transfer  additional  stresses.  The  stiffness 

258  recovery  can  be  approximated  with  the  parameter  A^,  defined  as  [26]: 


At  ^ 


VfEf 


VmEm  +  VfEf 


El  —  E2 

) - 

El 


(23) 


259  where  Vj  and  Vm  are  the  fiber  and  matrix  volume  fraction,  respectively.  Ef  and  E^a  are  the 

260  fiber  and  matrix  Young  modulus,  respectively,  b  is  and  adjustment  parameter  between  0,  if  the 

261  stiffness  is  completely  recovered,  and  1,  when  the  stiffness  recovery  is  only  due  the  matrix. 

262  Under  longitudinal  compressive  stresses  the  damage  variable  depends  of  the  damage  generated 

263  under  tension  and  compression  as: 


=  1  -  (1  -  di_)(l  - 


(24) 


264  The  damage  variable  (d^-)  is  assumed  to  be  approximated  by  an  exponential  law  as: 
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1 


(25) 


d 


* 

L- 


1  - 


rL- 


exp  [Al-  (1 


rL-)\ 


265  2.6.5  Shear  damage  law 

266  The  shear  stiffness  depends  on  the  longitudinal  and  transverse  damage.  It  is  assumed  that  if 

267  material  is  only  damaged  longitudinally  the  shear  damage  variable  takes  the  same  value  of  that 

268  of  the  longitudinal  tensile  variable  {ds  =  di+).  If  the  damage  is  due  the  transverse  damage,  the 

269  damage  variable  {d*g)  is  adjusted  with  the  transverse  critical  fracture  energy  in  mode  II.  The 

270  coupling  of  both  damage  variables  takes  the  following  form: 


ds  =  l-{l-d*s){l-dL+)  (26) 

271  Under  a  shear  test,  the  material  fails  in  transverse  direction.  Although  the  non- localized  damage 

272  is  important  in  many  hber  reinforced  plastics,  the  dehnition  of  a  constitutive  model  representing 

273  distributed  damage  and  plasticity  is  outside  the  scope  of  this  paper  and  will  be  addressed  in 

274  future  work.  Therefore,  an  exponential  law  that  enforce  softening  of  the  material  response  is 

275  proposed: 


d*s  =  l - exp  [As  (1  -  rr)] 

Tt 


276  Applying  the  crack  band  model,  the  parameter  As  results  in: 


As 


2iSl 

2Gi2Gs-iSl 


(27) 


(28) 


277  2.7  Tangent  eonstitutive  tensor  and  algorithm 


278  The  effective  computational  implementation  of  the  model  in  an  implicit  hnite  element  code 

279  requires  the  dehnition  of  the  tangent  stiffness  tensor  Ct: 


=  4  =  Cxe  (29) 

280  The  hrst  step  in  the  dehnition  of  the  tangent  stihness  tensor  is  the  calculation  of  the  time 

281  derivative  of  equation  (4): 


d=  Ce  with  C  =  H  ^  (I  -  M) 
282  where  the  matrix  M  is  dehned  as: 


(30) 
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283 


(7ii _ Odi  (7x1 _ dd\  (7x1 _ dd\  (7xi _ dd\  cr  xx _ dd-i 

{l-difEi  deii  {1-difEi  de22  {1-difEi  dess  {1-difEi  d'ns  (1-difEi  d'ns 

<722  dd2  722  dd2  722  dds  722  dds  722  dd2 

{l-dsfEs  deii  (i-dsfEs  de22  {l-dsfEs  dess  (l-dsfEs  d'ns  (l-dsfEs  dms 

M  =  733  dds  733  dds  733  dds  733  dds  733  dds 

{l-dsfE2  deii  (l-dsfEs  de22  {l-dsfEs  dess  {l-dsfEs  d'ns  (l-dsfEs  dms 

(7x2 _ Ode  (7x2 _ dde  <7x2 _ Ode  <7x2 _ OdQ  <7x2 _ OdQ 

il-defGi2  den  {l-defCis  dess  {l-ds)^Gis  dess  {l-dsfGi2  8712  (l-rf6)^Gi2  ^713 

(7X3 _ ddQ  (7X3 _ OdQ  (7X3 _ Od^  (7X3 _ OdQ  (7X3 _ Otie 

-  {^  —  dQ)^Gl2  (1  — cZ6)^Gi2  ^^22  (1  — d6)^Gi2  ^^33  {l  —  dQ)^Gl2  ^7l2  (1  — t/6)^Gx2  ^7l3 


Calculating  the  time  derivative  of  equations  (1)  and  (7)  and  using  equation  (30),  results  in; 


„C0.  ^CO 

^  =  T  ^  ^  Te  +  Te  (31) 

[  0  J  [ 0  Oj  [  0  0 

284  Consider  now  a  fixed  coordinate  system  that  coincides,  at  given  instant,  with  the  moving 

285  coordinate  system.  In  this  hxed  coordinate  system  the  transformation  matrix  is  the  identity 

286  T  =  I,  and  its  time  derivative  is:  T  =  ^  ^  d^  ■  j 

00000  0 
00000  1 
0  0  0  0  0  -1 
0  0  0  0  -1  0 
00010  0 
0  -2  2  0  0  0 

287  Using  (32)  in  (31),  the  tangent  stiffness  matrix  can  be  expressed  as  a  function  of  the  matrix  C: 


288  To  calculate  the  tangent  stiffness  tensor  (^  =  Ctc)  in  the  global  coordinate  system  it  is  necessary 

289  to  rotate  the  matrix  C  using  the  angle  6  as:  Ct  =  T^CT.  The  constitutive  model  results  in  an 

290  explicit  formulation  where  no  iterations  are  required  to  update  the  state  variables.  The  model 

291  developed  was  implemented  in  ABAQUS  non-linear  hnite  element  code  [46]  using  a  user-defined 

292  subroutine  UMAT,  according  to  the  following  algorithm; 
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1 

-  Strain  tensor  at  cnrrent  time  t 

e* 

2 

-  Main  transverse  direction 

9^ 

f  '\ 

3 

-  Main  transverse  strains 

r  l=TV 

loj 

4 

-  Ehective  stress  tensor 

d*  =  Ho 

5 

-  Loading  fnnctions 

(5-*) 

6 

-  Internal  variables 

/y«t  /  ^t-1  \ 

'  M  V  M  ^  ^M) 

7 

-  Damage  variables 

(^m) 

8 

-  Nominal  stress  tensor 

a*  = 

9 

-  Rotated  nominal  stress  tensor 

=  (T*)^  1 1 

lo  J 

10 

-  Tangent  constitntive  tensor 

C!^  =  (T*)^  C*T* 

294  3  Model  verification 

295  To  evaluate  the  accuracy  of  the  model  proposed,  the  numerical  predictions  are  compared  with 

296  published  experimental  results.  The  example  selected  for  the  assessment  of  the  accuracy  of 

297  the  model  corresponds  to  the  prediction  of  the  onset  and  accumulation  of  transverse  matrix 

298  cracks  and  of  hnal  failure  of  multidirectional  laminates.  This  is  a  relevant  scenario  because  the 

299  transverse  matrix  cracks  that  may  develop  at  low  strains  affect  the  thermomechanical  properties 

300  of  the  laminate,  create  paths  for  fnel  leakage,  and  may  trigger  other  failnre  mechanisms. 

301  Several  analytical  solntions  have  been  proposed  to  predict  the  effects  of  transverse  matrix 

302  cracks  in  the  thermomechanical  properties  of  mnltidirectional  laminates  [48]  -[55].  Generally, 

303  these  analytical  solntions  are  valid  for  simple  bonndary  conditions  and  applied  loads,  and 

304  for  sitnations  where  the  transverse  matrix  cracks  accnmnlate  in  a  central  90°  ply.  When  the 

305  cracked  90°  ply  is  placed  at  the  snrface  of  the  laminate  the  analytical  solntions  are  in  general 

306  no  longer  valid.  In  addition,  there  are  sitnations  where  delamination  develops,  either  preceding 

307  transverse  matrix  cracks  or  at  high  densities  of  such  cracks  [56].  The  model  proposed  here  is 

308  able  to  simulate  the  different  failnre  scenarios  where  the  analytical  solntions  are  no  longer  valid. 

309  3.1  Statistical  distribution  of  properties 

310  The  stress  held  of  an  nnnotched  laminate  is  nniform,  except  in  regions  close  to  the  free  edge 

311  where  a  three-dimensional  stress  held  occnrs.  This  means  that  the  damage  activation  fnnctions 

312  are  satished  in  more  than  one  point  simnltaneonsly  in  the  nnmerical  implementation  of  the 
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313  model.  To  overcome  this  difficulty,  and  to  enforce  the  localization  of  damage  in  a  plane,  a 

314  random  strength  held  is  created. 

315  The  inclusion  of  random  material  properties  in  the  hnite  element  mesh  is  an  active  research  topic 

316  [57]- [58].  The  dehnition  of  random  material  properties  must  take  into  account  the  characteristic 

317  length  of  the  hnite  element  and  should  be  applied  to  all  material  properties,  not  only  to  the  mean 

318  and  to  the  standard  deviation,  but  even  to  the  density  function  itself.  The  topic  of  statistical 

319  hnite  element  analysis  is  outside  the  scope  of  this  paper.  Therefore,  the  application  of  random 

320  material  properties  is  done  simply  by  dehning  a  random  strength  with  a  normal  distribution. 

321  A  simple  way  to  dehne  a  normal  distribution  is  by  means  of  the  Box-Muller  algorithm  [59]. 

322  The  Box-Muller  algorithm  states  that  if  a  and  b  are  uniformly  distributed  numbers  in  (0,1]  a 

323  standard  normally  distributed  random  variable  is  X  given  as:  X  =  V~2  In  a  cos(27r&) . 

324  If  the  maximum  and  minimum  strength  values,  and  respectively  are  known,  the 

325  mean  strength  value  is  /r  =  -1-  Y^^)/2  and  the  variance  is  given  as 

326  3(7  =  /i  —  Tj!?™  =  —  jj,.  The  random  normally  distributed  strength  variable  is:  Yt  = 

327  1/2  +  1/6  {y^^^  -  X. 

328  3.2  Kinematics  of  crack  growth 

329  The  predictions  of  the  onset  and  growth  of  transverse  matrix  cracks  use  a  [02/904]s  laminate, 

330  with  the  material  properties  shown  in  Tables  2  and  3.  A  typical  value  of  (7^+  for  glass  hbers  is 

331  used  in  the  simulations. 

332  [Table  2  about  here.] 

333  [Table  3  about  here.] 

334  Eight-node  solid  elements  (Abaqus  C3D8  elements)  are  used  in  the  fully  three-dimensional  nu- 

335  merical  models.  The  models  use  two  elements  through  the  thickness  of  each  ply.  The  specimen 

336  modeled  has  all  the  nodes  in  one  end  clamped,  whereas  the  other  extremity  is  subjected  to  an 

337  uniform  displacement. 

338  The  process  of  crack  propagation  in  a  [0„/90m]s  is  qualitatively  shown  in  the  Figures  3  and 

339  4.  The  damage  variable  d2  and  the  principal  strain  are  plotted  at  different  stages  of  cracking 

340  process  in  Figure  3.  As  expected,  damage  initiates  at  the  free  edge  of  the  laminate.  When 

341  increasing  the  external  load,  the  strain  localizes  and  transverse  matrix  cracks  develop.  Steps 

342  3-5  shown  in  Figure  3  represent  the  evolution  of  the  cracks  towards  the  center  of  the  laminate. 

343  The  microcracks  in  the  vicinity  of  through-the-thickness  cracks  unload  elastically.  The  process 

344  of  crack  grow  to  the  center  of  the  laminate  is  shown  in  Figure  4.  Due  the  conhning  effect  of 

345  outer  plies,  the  central  region  of  the  crack  advances  faster  than  the  region  in  the  vicinity  of  the 

346  adjoining  plies. 

347  [Fig.  3  about  here.] 

348  [Fig.  4  about  here.] 
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349  The  predictions  are  consistent  with  experimental  observations  [60]- [47]  that  show  that  matrix 

350  cracks  in  multiaxial  laminates  with  central  90°  plies  originate  at  the  free  edge  and  propagate 

351  through  the  specimen  width. 

352  When  [90„/0m]s  laminates  are  tested  in  tension,  the  matrix  cracks  appear  in  the  outer  plies 

353  in  an  antisymmetrical  distribution  [48].  Figure  5  shows  that  antisymmetrical  distribution  of 

354  cracks  is  correctly  captured  by  the  model. 

355  [Fig.  5  about  here.] 

356  Figure  6  represents  the  cracking  process  at  different  loading  stages.  The  presence  of  one  crack 

357  in  the  outer  ply  causes  the  loss  of  symmetry  of  the  laminate,  and  the  neutral  line  moves  away 

358  from  the  crack.  Therefore  the  laminate  bends  in  the  vicinity  of  the  crack.  The  strain  field  in 

359  the  location  where  a  crack  develops  is  shown  in  Figure  6. 

360  [Fig.  6  about  here.] 

361  The  previous  results  show  that  the  proposed  damage  model  is  able  to  qualitatively  represent 

362  the  process  of  matrix  cracking.  This  is  due  the  ability  of  the  hnite  element  model  to  describe 

363  the  kinematics  of  cracking  process  with  a  reasonable  accuracy. 

364  3.3  Effects  of  transverse  cracks  on  the  laminate  stiffness 

365  Varna  et  al.  [49]-[51]  presented  the  response  of  unnotched  {6  =  0°,  15°,  30°,  40°) 

366  glass/epoxy  laminates  loaded  in  tension.  The  reduction  of  the  laminate  Youngs  modulus  {Ex) 

367  and  Poisson  ratio  {vxy)  were  reported  as  a  function  of  the  laminate  strain  (Exx),  where  the 

368  direction  x  coincides  with  6  =  0°.  The  material  properties  are  obtained  from  [49]- [51]  and 

369  summarized  in  the  Tables  2  and  3. 

370  To  simulate  the  thermal  residual  strains  produced  due  the  curing  process,  a  temperature  change 

371  of  AT  =  —105  °C  is  applied.  The  ply  thickness  is  t  =  0.144mm.  The  other  parameters  required 

372  by  the  model  are  EIt-  =  0.4,  Xpo  =  600MPa,  Gf  =  Gf  =  30N/mm,  and  =  0. 

373  The  hnite  element  model  has  a  length  of  8mm  and,  taking  advantage  of  the  symmetry,  only 

374  the  one-half  of  the  laminate  thickness  is  modeled.  The  mesh  consists  of  C3D8  elements  with 

375  a  characteristic  length  of  0.072mm.  To  reduce  the  CPU  time  the  specimen  width  is  modeled 

376  with  only  one  row  of  hnite  elements.  Multi  point  constraints  are  introduced  in  the  model  as 

377  proposed  in  [65],  to  represent  a  state  of  generalized  plane  strain.  A  hrst  thermal  step  is  applied 

378  to  represent  the  curing  process.  A  prescribed  displacement  is  applied  to  one  end  of  the  specimen, 

379  while  the  opposite  end  is  clamped. 

380  No  transverse  damage  are  considered  in  the  outer  plies  when  6  ^  90°  plies  are  modeled.  In 

381  these  laminates,  matrix  cracks  develop  in  the  direction  described  by  the  ply  orientation,  and 

382  the  boundary  conditions  used  in  the  the  model  are  no  longer  valid. 

383  The  maximum  transverse  strength  used  is  =  200MPa.  Therefore,  the  mean  strength  is 

384  125MPa  and  the  standard  variation  is  25MPa.  The  strength  value  Sl  and  Yc  are  also  ran- 

385  domized  with  a  mean  value  of  180MPa  and  295MPa,  respectively.  The  standard  deviation  is 
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386  of  36MPa  and  59MPa,  respectively.  The  determination  of  the  random  strength  properties  was 

387  based  on  the  [02/904]^  laminate,  and  was  keep  constant  in  the  other  models.  The  inflnence  of  the 

388  random  properties  on  the  crack  pattern  is  small.  The  statistical  distribntion  mainly  inflnences 

389  the  onset  of  damage. 

390  It  should  be  noted  that  the  response  of  the  material  is  quite  sensitive  to  the  hnite  element  length. 

391  The  sensitivity  to  the  hnite  element  length  is  unrelated  to  the  common  problem  computing  the 

392  energy  dissipation  as  a  function  of  the  element  length;  this  problem  is  solved  by  virtue  of 

393  equation  16.  The  mesh  dependency  observed  is  due  to  the  relation  between  the  element  size 

394  and  the  statistical  distribution  used.  The  solution  of  this  problem  is  outside  the  scope  of  this 

395  work. 

396  For  the  [02/904]*  laminate,  the  predicted  failure  strain  is  =  0.022mm/mm  and  the  collapse 

397  of  the  laminate  occurs  when  hbre  fracture  develops  in  the  the  outer  0°  plies.  The  predicted  value 

398  is  smaller  than  the  one  obtained  using  classical  laminate  theory,  which  is  =  0.023mm/mm. 

399  This  fact  is  due  the  stress  concentration  that  the  matrix  cracks  cause  in  the  outer  plies,  which  in 

400  turn  leads  to  laminate  failure.  The  mean  stress-strain  response  up  to  laminate  failure  is  shown 

401  in  Figure  7. 

402  [Fig.  7  about  here.] 

403  Figures  8  to  11  show  the  internal  variable  tt  and  the  transverse  damage  variable  dT+  at  a  mean 

404  strain  =  0.01  mm/mm.  The  pattern  of  the  internal  variable  clearly  shows  the  localization  of 

405  the  deformations,  which  represent  transverse  cracks.  The  predicted  distribution  of  the  damage 

406  variable  indicates  that  delamination  develops  at  the  tip  of  the  transverse  matrix  crack. 

407  At  the  mean  strain  of  Sxx  =  O.Olmm/mm,  the  [02/004]*  laminate  shows  three  cracks,  cor- 

408  responding  to  a  crack  density  of  0.375  cracks/mm.  The  predicted  crack  density  of  the  other 

409  laminates  is  0.25  cracks/mm.  These  values  are  in  reasonable  agreement  with  the  values  reported 

410  by  Varna  [49],  which  are  0.34  cracks/mm  for  the  [O2/9O4]*  laminate,  0.28  cracks/mm  for  the 

411  [±15/904]*  laminate,  0.24  cracks/mm  for  the  [±30/904]*  laminate,  and  0.15  cracks/mm  for  the 

412  [±40/904]*  laminate. 

413  [Fig.  8  about  here.] 

414  [Fig.  9  about  here.] 

415  [Fig.  10  about  here.] 

416  [Fig.  11  about  here.] 

417  Figure  12  shows  the  longitudinal  and  transverse  tension  damage  variable  at  the  maximum  load 

418  and  after  laminate  failure.  It  is  possible  to  observe  that  laminate  failure  is  due  to  fiber  tensile 

419  fracture  of  the  adjoining  plies,  which  in  turn  is  triggered  by  the  transverse  matrix  cracks  that 

420  develop  in  the  inner  90°  layer. 

421  [Fig.  12  about  here.] 

422  The  reduction  in  the  laminate’s  Young  modulus  and  Poisson  ratio  are  shown  in  Figures  13  to  16. 

423  The  Young  modulus  is  calculated  using  the  predicted  load,  P,  the  predicted  end-displacement, 

424  (5,  the  specimen’s  cross-section  area.  A,  and  length,  L,  as  Ex  =  It  is  observed  that  the  crack 
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425  grows  in  the  transverse  direction  until  reach  a  critical  length  is  reached;  afterwards,  the  crack 

426  grows  unstably  with  an  associated  amount  of  spare  energy  that  produce  the  oscillations  that 

427  are  visible  in  the  response.  The  numerical  results  are  in  good  agreement  with  the  experimental 

428  results  reported  by  Varna  et  al.  [49]-[51]. 


429 


430 


431 


432 


[Fig.  13  about  here.] 
[Fig.  14  about  here.] 
[Fig.  15  about  here.] 
[Fig.  16  about  here.] 


433  3-4  Simulation  of  delamination 

434  The  model  is  further  validated  by  simulating  an  unidirectional  specimen  with  a  cut  across  the 

435  width  of  one  central  ply.  When  loaded  in  tension,  such  type  of  specimen  promotes  delamination 

436  growth  in  mode  II  [66]  before  the  outer  plies  fail  by  hber  tensile  fracture. 

437  The  specimen  under  investigation  was  tested  at  the  German  Aerospace  Centre  (DLR)  [67], 

438  and  it  consists  of  a  977-2  HTS  [O7]  CFRP  laminate  where  the  central  ply  was  cut  across  the 

439  entire  width  of  the  specimen.  The  relevant  material  properties  are  shown  in  Tables  4  and  5. 

440  The  nominal  ply  thickness  is  0.25mm. 

441  [Table  4  about  here.] 

442  [Table  5  about  here.] 

443  Five  specimens  were  tested  and  the  mean  value  of  the  remote  stress  at  delamination  propa- 

444  gation,  dehned  as  the  ratio  between  the  load  and  the  cross  section  area  of  the  specimen,  is 

445  1753MPa. 

446  The  modeling  strategy  used  here  consists  in  imposing  generalized  plane  strain  conditions  to 

447  Abaqus  C3D8  8-node  continuum  elements  by  applying  the  kinematic  relations  proposed  in 

448  [65].  The  nodes  at  one  end  of  the  specimen  are  clamped  and  a  displacement  in  the  longitu- 

449  dinal  direction  is  applied  to  the  nodes  at  the  other  end.  In  addition,  symmetry  along  the  1-2 

450  (longitudinal-transverse)  midplane  of  the  specimen  is  imposed.  The  ply  cut  was  simulated  by 

451  a  line  of  elements  with  all  the  damage  variables  set  to  one.  The  mesh  of  the  specimen  under 

452  investigation  is  shown  in  Figure  17. 

453  [Fig.  17  about  here.] 

454  Figure  18  shows  the  predicted  relation  between  the  end  displacement  of  the  specimen  and  the 

455  remote  stress.  The  propagation  of  delamination  corresponds  to  the  horizontal  plateau  predicted 

456  by  the  model. 

457  [Fig.  18  about  here.] 

458  Figures  19  and  20  show  the  the  sequence  of  the  failure  mechanisms  that  occur  in  the  specimen. 


18 


459 


[Fig.  19  about  here.] 
[Fig.  20  about  here.] 


460 

461  It  is  observed  that  the  ply  cut  triggers  a  delamination  between  the  central  ply  and  the  adjoining 

462  plies,  that  propagates  in  mode  II  along  the  length  of  the  specimen.  The  specimen  is  able  to 

463  sustain  increasing  loads  until  it  completely  fails  as  a  result  of  the  fiber  fracture  in  the  adjoining 

464  plies.  This  sequence  of  events  was  also  observed  in  the  experiments.  In  addition,  the  predicted 

465  remote  stress  at  delamination  propagation,  1782. 4MPa,  correlates  well  with  the  mean  value 

466  measured  in  the  experiments,  1753MPa. 


467  4  Conclusions 

468  A  fully  three-dimensional  continuum  damage  model  able  to  predict  the  different  failure  mech- 

469  anisms  that  may  occur  in  laminated  composites  was  proposed.  The  constitutive  model  is  for- 

470  mulated  in  the  formalism  of  the  thermodynamics  of  irreversible  processes  and  its  numerical 

471  implementation  ensures  a  mesh-independent  prediction  of  energy  dissipation  by  using  the  crack 

472  band  model. 

473  The  preliminary  validation  examples  indicate  that  the  model  is  able  to  capture  the  kinematics  of 

474  the  propagation  of  transverse  matrix  cracks  for  quasi-isotropic  laminates  with  general  locations 

475  of  the  90°  plies.  The  comparison  between  the  model  predictions  and  published  experimental 

476  data  show  that  the  model  is  able  to  accurately  predict  the  relation  between  the  applied  strain 

477  and  the  residual  stiffness  of  quasi-isotropic  laminates,  hnal  failure  of  the  laminates,  as  well  as 

478  the  effect  of  the  stiffness  of  the  adjoining  sublaminates  on  the  density  of  transverse  matrix 

479  cracks.  In  addition,  the  model  is  able  to  represent  at  the  constitutive  level  both  delamination 

480  and  transverse  matrix  cracks  as  well  as  the  interaction  between  these  two  failure  mechanisms 

481  that  occurs  at  in  glass-epoxy  quasi-isotropic  laminates  at  high  applied  strains. 

482  Based  on  the  results  of  the  simulation  of  a  CFRP  specimen  with  a  central  cut  ply  it  is  concluded 

483  that  the  model  accurately  represents  delamination  onset  and  propagation  in  mode  II,  and  the 

484  hnal  fracture  of  the  laminate  as  a  result  of  hber  fracture.  It  should  be  emphasized  that  this 

485  sequence  of  events  is  predicted  without  recurring  to  special  purpose  cohesive  elements,  but 

486  using  an  appropriate  constitutive  model  for  the  composite  ply. 

487  The  future  research  of  the  authors  will  include  a  further  validation  of  the  model  presented  in  this 

488  paper  for  loading  scenarios  that  trigger  compression  and  shear-dominated  failure  mechanisms. 
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Fig.  1.  Transverse  damage  activation  function,  Ft  =  0. 
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Fig. 


.  Propagation  of  transverse  crack. 
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Fig.  5.  Distribution  of  transverse  matrix  cracks  in  a  [90/0]^  laminate.  Deformed  scale:  30x. 
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Fig.  6.  Cracking  process  in  a  [90/0]  s  laminate.  Deformed  scale:  30x. 
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Fig.  8.  Internal  variable  and  transverse  tension  damage  variable  {dj'+)  for  a  [02/904]^  laminate. 
Mean  laminate  deformation  of  exx=0-0^-  Deformed  scale:  lOx. 
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Fig.  9.  Internal  variable  rrp  and  transverse  tension  damage  variable  dT+  for  a  [±15/904]^  laminate. 
Mean  laminate  deformation  of  exx=0-0^-  Deformed  scale:  lOx. 
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Fig.  10.  Internal  variable  and  transverse  tension  damage  variable  dT+  for  a  [±30/904]^  laminate. 
Mean  laminate  deformation  of  exx=0-01.  Deformed  scale:  lOx. 
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Fig.  11.  Internal  variable  and  transverse  tension  damage  variable  dT+  for  a  [±40/904]^  laminate. 
Mean  laminate  deformation  of  exx=0-01.  Deformed  scale:  lOx. 
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Fig.  12.  Longitudinal  di-^.  and  transverse  tension  damage  variable  dT+  for  a  [02/904]^  laminate  at 
failure  load  and  when  the  specimen  is  completely  broken.  Deformed  scale  lOx. 
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Cflass/Epoxy  laminate 


Fig.  13.  Ex  and  Vxy  as  functions  of  the  applied  strain  for  a  [O2/9O4]  laminate.  Experimental  results 
from  Varna  et  al.  [49]-[51]. 
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[±  15/90^]^  Glass/Epoxy  laminate 


Fig.  14.  Ex  and  Uxy  as  functions  of  the  applied  strain  for  a  [±15/904]s  laminate.  Experimental  results 
from  Varna  et  al.  [49]-[51]. 
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[±  30/90^]^  Glass/Epoxy  laminate 


Fig.  15.  Ex  and  Uxy  as  functions  of  the  applied  strain  for  a  [±30/904]^  laminate.  Experimental  results 
from  Varna  et  al.  [49]-[51]. 
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[±  40/90^]^  Glass/Epoxy  laminate 


Fig.  16.  Ex  and  Vxy  as  functions  of  the  applied  strain  for  a  [±40/904]^  laminate.  Experimental  results 
from  Varna  et  al.  [49]-[51]. 
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Fig.  17.  Mesh  of  the  unidirectional  specimen. 
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Remote  Stress  (MPa) 


2000-1 


Fig.  18.  Predicted  relation  between  the  remote  stress  and  the  end  displacement. 
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Fig.  19.  Damage  variable  dg  at  1782. 4MPa.  Deformed  scale:  2x. 
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Fig.  20.  Damage  variable  di  at  peak  load,  1829MPa.  Deformed  scale:  2x. 
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Table  1 

List  of  symbols. 


Cij 

Strain  in  the  material  coordinate  system. 

^ij 

Nominal  stress  in  the  material  coordinate  system. 

£ij 

Principal  strain  in  the  transversely  isotropic  plane. 

O'ij 

Principal  nominal  stress  in  the  transversely  isotropic  plane. 

O'ij 

Principal  effective  stress  in  the  transversely  isotropic  plane. 

T{e) 

Transformation  matrix  between  and  £ij. 

i’ 

Complementary  free  energy  per  unit  volume. 

AT 

Difference  in  temperature. 

AM 

Difference  in  moisture  content. 

a 

Coefficients  of  thermal  expansion. 

P 

Coefficients  of  hygroscopic  expansion. 

di 

Active  damage  variables  (I  =  1,  2,  3,  6). 

cLm 

Damage  variables  (M  =  L+,  L—,T+,T—,  S). 

El^  E2,  (712,  J^12,  2^23 

Ply  engineering  elastic  constants. 

Ex ,  Vxy 

Laminate  engineering  elastic  constants. 

H 

Compliance  tensor. 

Yi 

Thermodynamic  forces  associated  to  active  damage  variable. 

Rate  of  dissipation. 

En 

Damage  activation  functions  {N  =  L+,L—,T). 

4>n 

Loading  functions  {N  =  L+,  L—,T). 

tn 

Elastic  domain  thresholds  {N  =  L+,  L—,  T). 

Xt,  Xc,  Tt,  Tc,  Sl 

Uniaxial  strengths. 

1 

Characteristic  element  size. 

Gm 

Fracture  toughness  (M  =  L+,  L—,T+,  S). 

9m 

Dissipated  energy  per  unit  volume  (M  =  L+,  L—,T+,  S). 

Hj'— 

Transverse  compression  damage  parameter. 

Xpo 

Pull-out  strength. 

r-<E 

^L+ 

Pull-out  fracture  energy. 

/1± 

Longitudinal  tension-compression  coupling  parameter. 
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Table  2 

Elastic  properties  of  glass-epoxy  [49]- [51]. 


El  (MPa) 

E2  (MPa) 

Gi2  (MPa) 

1^12  1^23  an  (1/°C)  022  (1/°C) 

44730 

12760 

5800 

0.28  0.42  8.6x10“^^  22.1x10“® 
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Table  3 

Strength  and  fracture  energy  of  glass-epoxy. 


Xt  (MPa) 

Yt  (MPa) 

Sl  (MPa) 

Gt+  (N/mm) 

Gs  (N/mm) 

Gl+  (N/mm) 

1060 

50 

72 

0.4 

0.8 

60 
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Table  4 

Elastic  properties  of  977-2  HTS. 


El  (MPa) 

E2  (MPa) 

Gi2  (MPa)  1^12  1^23 

144000 

7500 

5030  0.29  0.50 
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Table  5 


Strengths  and  fracture  toughness  of 

977-2  HTS. 

Xt  (MPa) 

Yt  (MPa) 

Sl  (MPa) 

Gs  (N/mm) 

2290 

47 

67 

1.55 
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